Pricing II: Aplicación de R a la tarificación

Author

Carolina Armella, Ágatha del Olmo y María Yu García

Published

October 18, 2024

1. Introducción

El propósito de este análisis es predecir el número esperado de reclamaciones que realizará una póliza a partir de la experiencia de siniestralidad histórica de una cartera de seguros. Este estudio se basa en la implementación de modelos lineales generalizados (GLM), que permiten realizar predicciones precisas manteniendo un alto grado de interpretabilidad.

Para este caso, utilizamos un conjunto de datos reales proporcionado por una compañía de seguros, que recoge información sobre 30,000 pólizas. En el presente análisis, para simplificar y reducir el tiempo de procesamiento, trabajamos con una muestra aleatoria de la base de datos original.

El análisis incluye una exploración detallada de los datos (Exploratory Data Analysis o EDA) y la implementación de modelos de aprendizaje automático para evaluar diferentes metodologías. Además, dado el contexto regulatorio de la industria de seguros, donde la transparencia y la explicabilidad son fundamentales, se priorizan los modelos que cumplen con estos requisitos, por eso no vamos a usar modelos de tipo caja negra.

2. Metodología

Para implementar un sistema automatizado de evaluación crediticia (credit scoring), es fundamental contar con modelos que permitan identificar con precisión los riesgos asociados a los préstamos. Además, la normativa bancaria en muchos países, incluido España, exige que las instituciones expliquen a los solicitantes los motivos para aceptar o rechazar una solicitud de crédito. Por ello, los modelos empleados deben ser transparentes. Herramientas más avanzadas como las redes neuronales o los modelos de máquinas de soporte vectorial (SVM) son consideradas “cajas negras” y no son fáciles (y a veces son imposibles) de interpretar y por lo tanto no cumplen con este requisito de transparencia.

2.1 Exploración de los datos

La base de datos corresponde a una extracción aleatoria de la BBDD de experiencia de siniestralidad de una compañía de seguros y consta de 19 variables capturando distintas características de las pólizas.

Reducción del dataset

Hemos decidido reducir la base de datos a 1000 observaciones en vez de las 27378 originales para agilizar la convergencia por problemas de rendimiento.

datos <- read.csv("D:\\3 curso\\TP\\Grupo 2.csv", sep=";")

datos <- datos[1:1000,] # reducimos el training set para que tarde menos

dim(datos)
[1] 1000   19
str(datos)
'data.frame':   1000 obs. of  19 variables:
 $ npoliza     : int  417952 522479 911639 587785 1117004 560994 414660 1095007 1074745 171239 ...
 $ iniciop     : int  19990101 20010712 20000101 20011110 20011205 20010923 20011218 20010818 20010101 20000101 ...
 $ finalp      : int  20000101 20020101 20000120 20020101 20020101 20020101 20020101 20020101 20010530 20001204 ...
 $ potencia    : int  136 57 75 68 115 100 96 75 60 90 ...
 $ valorvehi   : int  4950 1900 2000 1500 2650 2450 3350 2000 1800 1578 ...
 $ antivehi    : int  2 12 1 14 6 4 9 1 2 11 ...
 $ tipovehi    : chr  "100" "100" "100" "100" ...
 $ plazasvehi  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ edad        : int  40 49 41 40 47 37 42 47 49 31 ...
 $ anticarnet  : int  12 31 20 20 24 17 23 25 10 12 ...
 $ sexocondu   : chr  "V" "V" "V" "V" ...
 $ anticia     : int  2 4 0 4 1 4 5 1 0 8 ...
 $ codigopostal: int  28003 20500 7702 37100 8291 28945 8006 38380 43840 35110 ...
 $ usovehi     : int  11000 11100 11000 11000 11000 11000 11000 11000 11000 11000 ...
 $ cuantiasi   : int  223849 0 0 0 0 0 0 0 0 111116 ...
 $ numerosi    : int  2 0 0 0 0 0 0 0 0 2 ...
 $ siniestro   : int  1 0 0 0 0 0 0 0 0 1 ...
 $ prov        : int  28 20 7 37 8 28 8 38 43 35 ...
 $ ccaa        : int  15 1 7 5 8 15 8 17 8 17 ...

En este caso, nuestra variable respuesta es numerosi, que indica el número de siniestros reportados por cada póliza. Esta variable es crucial porque el objetivo del modelo es predecir el número esperado de reclamaciones en función de las características de la póliza y del asegurado.

La variable nploiza es simplemente un identificador único de cada póliza, y no contiene información relevante o explicativa para la predicción de siniestros, por lo tanto la podemos eliminar.

Además, hay dos variables que nos dan información adicional. Podemos sustituir iniciop y finalp por la duración de la póliza en días o meses. Este predictor podría ser relevante teniendo en cuenta la exposición al riesgo. En cuanto al código postal, aunque en apariencia es solo una etiqueta numérica, contiene información geográfica valiosa, y para hacer más fácilmente interpretable este predictor extraemos la comunidad autónoma ya que no contamos con los metadatos de ccaa para saber a qué referencia cada número, por lo tanto esa variable la eliminamos.

2.2 Preprocesamiento de los datos

Eliminamos por las razones indicadas la variable npoliza.

# Descartamos el número de póliza
datos <- subset(datos, select = -c(npoliza))
Mostrar código
datos <- datos %>%
  mutate(region = substr(as.character(codigopostal), 1, 2)) %>%
  mutate(comunidad = case_when(
    region %in% c("01", "48", "20") ~ "País Vasco",
    region %in%  c("31") ~ "Navarra",
    region %in% c("03", "12", "46") ~ "Comunidad Valenciana",
    region %in% c("21", "41", "14", "29", "11", "18", "23", "04") ~ "Andalucía",
    region %in% c("08", "17", "25", "43") ~ "Cataluña",
    region %in% c("15", "27", "32", "36") ~ "Galicia",
    region %in% c("33") ~ "Asturias",
    region %in% c("39") ~ "Cantabria",
    region %in% c("09", "42", "47", "34", "37", "05", "24", "40", "49") ~ "Castilla y León",
    region %in% c("45", "13", "16", "19", "02") ~ "Castilla-La Mancha",
    region %in% c("22", "44", "50") ~ "Aragón",
    region %in% c("30") ~ "Murcia",
    region %in% c("28") ~ "Madrid",
    region %in% c("07") ~ "Islas Baleares",
    region %in% c("35", "38") ~ "Canarias",
    region %in% c("06", "10") ~ "Extremadura",
    region %in% c("26") ~ "La Rioja",
    region %in% c("51", "52") ~ "Ceuta y Melilla"
  ))

datos <- datos[, -19]

Comprobamos que no haya NULL en la variable creada ya que hemos contemplado todas las posibles combinaciones de códigos postales y sus comunidades autónomas según el INE en la fórmula.

sum(is.na(datos$comunidad))
[1] 260

Como vemos al crear la nueva variable comunidad autónoma, nos encontramos con muchos NA, por lo tanto, o la base de datos es completamente inventada o no sigue el sistema de códigos postales oficial, lo cual hace imposible su uso. Por lo tanto, eliminamos la variable creada.

datos <- datos[, -19]

Muchas de las variables no tienen la clase que les correspondería. Primero las variables iniciop y finalp las pasamos a formato fecha. También creamos la variable duración póliza con la diferencia entre estas dos, que nos indica cuánto ha durado la póliza en días.

# Convertimos en formato fecha: as.Date(x, "%m/%d/%Y")  
datos$iniciop <- paste0(substr(datos$iniciop, 5, 6), "/",
                        substr(datos$iniciop, 7, 8), "/",
                        substr(datos$iniciop, 1, 4) )

datos$iniciop <- as.Date(datos$iniciop , "%m/%d/%Y")

datos$finalp <- paste0(substr(datos$ finalp, 5, 6), "/",
                       substr(datos$ finalp, 7, 8), "/",
                       substr(datos$ finalp, 1, 4) )

datos$finalp <- as.Date(datos$ finalp , "%m/%d/%Y")

datos$duracion_poliza <- as.numeric(datos$finalp - datos$iniciop)

A veces también es necesario convertir variables a tipo factor. Los factores se usan comúnmente para representar variables categóricas, es decir, variables que toman un número limitado de valores. En nuestro caso:

  • tipovehi: Representar diferentes tipos de vehículos, en concreto 7.
  • codigopostal: Se convierte en un factor, ya que representa códigos postales, que son categorías y no valores numéricos continuos.
  • prov: Similar, representa provincias que se pueden tratar como categorías.
  • ccaa: Representa la comunidad autónoma (en España), lo que también es una variable categórica.
datos$tipovehi <- as.factor(datos$tipovehi)

datos$codigopostal <- as.factor(datos$codigopostal)

datos$prov <- as.factor(datos$prov)

datos$ccaa <- as.factor(datos$ccaa)

str(datos)
'data.frame':   1000 obs. of  19 variables:
 $ iniciop        : Date, format: "1999-01-01" "2001-07-12" ...
 $ finalp         : Date, format: "2000-01-01" "2002-01-01" ...
 $ potencia       : int  136 57 75 68 115 100 96 75 60 90 ...
 $ valorvehi      : int  4950 1900 2000 1500 2650 2450 3350 2000 1800 1578 ...
 $ antivehi       : int  2 12 1 14 6 4 9 1 2 11 ...
 $ tipovehi       : Factor w/ 7 levels "100","120","150",..: 1 1 1 1 3 1 2 1 1 1 ...
 $ plazasvehi     : int  5 5 5 5 5 5 5 5 5 5 ...
 $ edad           : int  40 49 41 40 47 37 42 47 49 31 ...
 $ anticarnet     : int  12 31 20 20 24 17 23 25 10 12 ...
 $ sexocondu      : chr  "V" "V" "V" "V" ...
 $ anticia        : int  2 4 0 4 1 4 5 1 0 8 ...
 $ codigopostal   : Factor w/ 747 levels "1000","1003",..: 439 371 101 587 167 462 112 613 666 552 ...
 $ usovehi        : int  11000 11100 11000 11000 11000 11000 11000 11000 11000 11000 ...
 $ cuantiasi      : int  223849 0 0 0 0 0 0 0 0 111116 ...
 $ numerosi       : int  2 0 0 0 0 0 0 0 0 2 ...
 $ siniestro      : int  1 0 0 0 0 0 0 0 0 1 ...
 $ prov           : Factor w/ 50 levels "1","2","3","4",..: 28 20 7 37 8 28 8 38 43 35 ...
 $ ccaa           : Factor w/ 17 levels "1","2","3","4",..: 15 1 7 5 8 15 8 17 8 17 ...
 $ duracion_poliza: num  365 173 19 52 27 100 14 136 149 338 ...

Echamos un vistazo a las tablas de las variables no continuas para ver su frecuencia:

Creamos ahora una nueva variable llamada usovehi2 que clasifica la variable numérica usovehi en dos grupos, los que tienen un valor menor a 2000 que se pasan a llamar G1 y los que tienen un valor mayor, que se llamarán G2. Después eliminamos la variable original.

datos$usovehi2 <- ifelse(datos$usovehi < 20000, "G1", "G2")

# Descartamos uso de vehiculo, poca varianza
datos <- subset(datos, select = -c(usovehi))

Observamos la matriz de correlaciones, que nos indica el coeficiente de correlación entre todos los pares de variables numéricas de la base de datos.

Coeficiente de correlación

El coeficiente de correlación mide la intensidad y dirección de la relación lineal entre dos variables:

  • \(r=1\): Correlación positiva perfecta, a medida que una variable aumenta, la otra también aumenta proporcionalmente.
  • \(r=−1\): Correlación negativa perfecta, a medida que una variable aumenta, la otra disminuye proporcionalmente.
  • \(r=0\): No hay correlación lineal, las variables no tienen una relación lineal directamente.
numericas <- datos[sapply(datos, is.numeric)]

# Calcular la matriz de correlación
correlaciones <- cor(numericas, use = "pairwise.complete.obs")

# Crear un gráfico de correlación
corrplot(correlaciones, method = "circle", type = "upper", 
         tl.col = "black", tl.srt = 45, col = colores)

Como podemos ver, lógicamente las variables están perfectamente correlacionadas de forma directa con ellas mismas. Pero podemos ver otras relaciones interesantes, para empezar hay una correlación bastante alta entre la variable objetivo numerosi y las variables cuantiasi y siniestro. Descartamos siniestro y cuantiasi como predictores.

cor(datos$siniestro, datos$numerosi)
[1] 0.7607234
cor(datos$cuantiasi, datos$numerosi)
[1] 0.3509457
datos <- subset(datos, select = -c(cuantiasi, siniestro))

Esto se debe a que la variable siniestro es un resultado de un evento ya ocurrido, y utilizarla para predecir el número de siniestros futuros es como usar la variable respuesta en el modelo de entrenamiento, lo que resulta en un modelo irrealista o sobreajustado por redundancias. Esto lo acabamos de ver reflejado en la matriz de correlaciones (\(r=0.76\)).

Además, al estar cuantiasi también altamente correlacionada con la variable de respuesta (\(r=0.35\)), usarla como predictor puede causar multicolinealidad por redundancia, lo que afecta negativamente la interpretación y precisión de los coeficientes del modelo.

El gráfico de correlaciones también nos ha mostrado una muy alta correlación entre potencia y valorvehi:

cor(datos$valorvehi, datos$potencia)
[1] 0.8207117

Esto nos lleva a un alto riesgo de multicolinealidad, que ocurre cuando dos o más variables predictoras están altamente correlacionadas entre sí, lo que puede causar problemas en los modelos estadísticos, especialmente en los modelos de regresión, ya que puede dificultar la interpretación de los coeficientes tal como hemos mencionado con cuantiasi. Por lo tanto, eliminamos la variable de menos interés interpretativo para el estudio, valorvehi, ya que la potencia es algo más objetivo que no depende de la economía.

datos <- subset(datos, select = -valorvehi)

Por último, el gráfico ha indicado una alta correlación entre edad y anticarnet.

cor(datos$anticarnet, datos$edad)
[1] 0.8174436

Si el estudio se enfocara en la relación entre estas dos variables para analizar el efecto de sacarse el carnet antes o después sería interesante mantener las dos, pero sabiendo que en la inmensa mayoría de casos las personas se sacan el carnet alrededor de los 18 años podríamos inferir de forma bastante acertada a partir de la edad la antigüedad del carnet. Por lo tanto, eliminamos la antigüedad de carnet.

datos <- subset(datos, select = -anticarnet)

Ahora podemos volver a ver el gráfico de correlaciones para comprobar que no quedan altas correlaciones:

numericas <- datos[sapply(datos, is.numeric)]

correlaciones <- cor(numericas, use = "pairwise.complete.obs")

corrplot(correlaciones, method = "circle", type = "upper", 
         tl.col = "black", tl.srt = 45, col = colores)

Al reducir el número de variables considerando las correlaciones, obtenemos una base de datos más manejable e interpretable, lo que facilita un análisis más eficiente y una convergencia más rápida en los modelos.

2.3 Separación de los datos en Training set y Test set

Partimos al azar base en dos: 66.7% entrenamiento (667 casos) y 33.33% test (333 casos).

Separación del dataset
  • Training set: Este conjunto se utiliza para entrenar el modelo. Contiene ejemplos históricos con características de préstamos y la información de si fueron pagados o no. El modelo analiza estos datos para identificar patrones que permitan clasificar a futuros solicitantes como buenos pagadores o posibles impagados.

  • Test set: Este segundo conjunto se emplea para evaluar la efectividad del modelo en datos que no ha visto antes. Sirve para comprobar si el modelo es capaz de reconocer patrones en casos nuevos y, por tanto, si sería útil para predecir con precisión en situaciones reales.

Para separar la base de datos se establece una semilla para garantizar que los resultados sean reproducibles. Luego, los datos se ordenan pseudo-aleatoriamente, se calcula el número de muestras para el conjunto de entrenamiento y de test, y se seleccionan los datos para ambos de la database original.

Después, se genera una tabla de frecuencias que muestra la distribución del número de reclamaciones en el conjunto de entrenamiento.

# Permutamos los datos al azar
set.seed(216514)
azar <- sample(1:nrow(datos))
casos <- round(nrow(datos)*0.667)

train <- datos[azar[1:casos],]
test <- datos[azar[(casos+1):nrow(datos)],]

Tabla.Training <- table(train$numerosi)
Tabla.Test <- table(test$numerosi)
Tabla.Training %>%
  kbl(col.names = c("Número Reclamaciones", "Frecuencia"),
      align = rep('c', 2)) %>%
  kable_styling()
Número Reclamaciones Frecuencia
0 489
1 115
2 35
3 11
4 9
5 5
6 2
8 1

Realizamos la misma tabla para el test set:

Tabla.Test %>%
  kbl(col.names = c("Número Reclamaciones", "Frecuencia"),
      align = rep('c', 2)) %>%
  kable_styling()
Número Reclamaciones Frecuencia
0 244
1 66
2 14
3 5
4 1
5 1
6 1
8 1

En ambos casos la tabla sugiere que la mayoría de las pólizas en el conjunto de datos tienen 0 o 1 reclamación. Las otras categorías tienen menos frecuencia, lo que indica que las reclamaciones mayores son poco frecuentes en esta muestra, lo cual era de esperar, e interesa mantener a la empresa.

La técnica de separación usada es una de las más comunes y simples, pero puede no ser la mejor opción en todos los casos, especialmente si la distribución de las clases de los datos están desbalanceadas como acabamos de comprobar. Teniendo en cuenta esta característica de nuestra database, podemos usar el muestreo estratificado. Este método asegura que la distribución de las clases en los conjuntos de entrenamiento y test sea representativa de la distribución en el conjunto original. Es útil para evitar sesgos en la evaluación del modelo, especialmente cuando las clases menos frecuentes son importantes (que en nuestro caso lo son ya que son las que más afectan negativamente en la economía de la empresa).

3. Análisis Exploratorio de Datos

3.1 Variable respuesta

Visualizamos la frecuencia de número de reclamaciones a lo largo del tiempo.

Mostrar código
barplot(table(train$numerosi),
        col = colores,                    
        border = "white",                
        main = "Frecuencia de Número de Reclamaciones", 
        xlab = "Número de Reclamaciones", 
        ylab = "Frecuencia", 
        las = 1,                          
        cex.names = 0.9,                 
        cex.axis = 0.8,                  
        cex.main = 1.2)             

Una de las propiedades más importantes de la distribución de Poisson es que su media es igual a su varianza. Ambos están representados por el parámetro λ (lambda). Esta propiedad es única y ayuda a identificar si un conjunto de datos sigue una distribución de Poisson. Esta igualdad implica que, a medida que aumenta el número esperado de sucesos, también lo hace la variabilidad en el número real de sucesos.

Por lo tanto, podemos encontraar sobredispersión si la varianza de la variable de respuesta es significativamente mayor que su media, lo cual es incompatible con los supuestos del modelo Poisson como acabamos de ver.

mean(train$numerosi)
[1] 0.4482759
var(train$numerosi)
[1] 0.9383867

La media es 0.45, mientras que la varianza es 0.94. En un modelo Poisson ideal, la media y la varianza deberían ser iguales.

La varianza es notablemente mayor que la media, lo que indica sobredispersión. Esto significa que los datos probablemente no sean bien modelados por un modelo Poisson estándar, ya que este asume que la media y la varianza son iguales. Esta igualdad implica que, a medida que aumenta el número esperado de sucesos, también lo hace la variabilidad en el número real de sucesos.

El exceso de ceros ocurre cuando la proporción de ceros en la variable de respuesta es mucho mayor de lo esperado bajo una distribución Poisson.

table(train$numerosi)[1] / nrow(train)
        0 
0.7331334 

Más del 73% de las observaciones tienen un valor de cero en número de reclamaciones, lo cual esperaríamos porque indica que muchos pólizas no tienen reclamaciones, pero en términos técnicos nos puede llevar a complicaciones. Este exceso de ceros no es consistente con un modelo Poisson estándar, que generaría ceros con menor frecuencia dado el valor esperado de 0.39.

Concluimos que el modelo Poisson estándar no es adecuado. Necesitarás un modelo que pueda manejar sobredispersión y exceso de ceros como por ejemplo los modelos de inflación cero. Estos utilizan una técnica estadística adaptada a los datos de recuento que tienen un exceso de resultados cero, lo cual es ideal para nuestro caso.

Ahora, representamps el número de reclamaciones ajustadas según un modelo de Poisson, el número de reclamaciones que el modelo de Poisson predice que deberíamos tener para cada número de reclamaciones frente al número de reclamaciones observadas en los datos reales.

Mostrar código
predict0 <- dpois(0:10, mean(train$numerosi)) * nrow(train)
observado <- table(factor(train$numerosi, levels = 0:10))

barplot(rbind(observado, predict0), las = 1,
        beside = TRUE, names.arg = c(0:9, ">=10"),
        col = c(colores[1], colores[5]),
        border = "white",                
        main = "Distribución Poisson", 
        xlab = "Número de Reclamaciones", 
        ylab = "Frecuencia", 
        las = 1,                          
        cex.names = 0.9,                 
        cex.axis = 0.8,                  
        cex.main = 1.2)

Como podemos ver, el modelo Poisson predice que deberíamos tener menos pólizas con 0 reclamaciones que las observadas, pero más de 1 reclamación, lo cual implicaría una variabilidad más constante con la media.

Comparamos con una Binomial Negativa:

Mostrar código
fit.bn <- goodfit(train$numerosi, type = "nbinomial", 
                  method = "MinChisq")

barplot(rbind(observado, fit.bn$fitted[1:11]), las = 1,
        beside = TRUE, names.arg = c(0:9, ">=10"),
        col=c(colores[1], colores[5]),
        border = "white",                
        main = "Distribución Binomial Negativa", 
        xlab = "Número de Reclamaciones", 
        ylab = "Frecuencia", 
        las = 1,                          
        cex.names = 0.9,                 
        cex.axis = 0.8,                  
        cex.main = 1.2)                

Claramente la distribución Binomial Negativa se ajusta de manera más adecuada a nuestra base de datos, ya que las frecuencias ajustadas se alinean mucho mejor con las observadas, reflejando de manera más precisa la estructura de los datos.

3.2 Relaciones variable respuesta-predictores

3.2.1 Tiempo de exposición al riesgo

El tiempo de exposición al riesgo se debe tratar de la siguiente manera según el contexto:

  • Como un offset: es el enfoque más común en modelos de conteos (Poisson, Binomial Negativa).

  • Como predictor: este enfoque es menos común, pero útil si creemos que el efecto del tiempo varía entre los casos.

  • Como peso: si queremos dar más importancia a los casos con más tiempo de exposición al riesgo.

Por lo tanto, para modelos de eventos o conteos, el tratamiento más adecuado es generalmente como un offset.

# Tiempo de exposición (Predictor, Peso u Offset)
summary(train$duracion_poliza)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    79.0   179.0   178.4   273.0   366.0 

Como el tiempo de exposición afecta la tasa de siniestros, modelarlo como un offset nos permitiría ajustar por esta variable sin estimar un coeficiente independiente para ella. Al tratarlo como offset, estaríamos ajustando el modelo para que la tasa de eventos sea proporcional al tiempo de exposición, pero sin incluirla como un predictor adicional cuyo coeficiente se estimaría.

Es posible que algunos contratos tengan una duración extremadamente corta, de menos de un día, como un seguro de alquiler que se contrata y se extingue el mismo día. En este caso, la “exposición al riesgo” podría ser efectivamente cero, ya que no hay un período suficiente de tiempo para que ocurra un evento asegurado (como un accidente o siniestro).

sum(datos$duracion_poliza==0)
[1] 24
sum(train$duracion_poliza==0)
[1] 13
sum(test$duracion_poliza==0)
[1] 11

En nuestro caso, el conjunto de entrenamiento contiene 13 cero días al azar de los 24 registros con ceros originalmente. De la misma forma, en el conjunto de prueba hay 11 casos con duración de póliza igual a cero.

En ambos conjuntos hay registros con duración de póliza de 0 días, lo cual podría ser problemático, así que los convertimos en 1, y lo repetimos para la posibilidad de que haya valores erróneos menores a 0 que se deberían a un error de imputación de los datos originales.

test$duracion_poliza[test$duracion_poliza == 0] <- 1
test$duracion_poliza[test$duracion_poliza < 0] <- 1
train$duracion_poliza[train$duracion_poliza == 0] <- 1
train$duracion_poliza[train$duracion_poliza < 0] <- 1

Para entender mejor la distribución de la duración de las pólizas en nuestro conjunto de datos, hemos realizado un gráfico de densidad. Este gráfico nos permite observar cómo se distribuyen los valores de duración de las pólizas, destacando posibles tendencias o patrones en la variable.

Mostrar código
ggplot(train, aes(x = duracion_poliza)) +
  geom_density(fill = colores[1], color = colores[5], alpha = 0.6) + 
  theme_minimal() +
  labs(title = "Distribución de la Densidad de la Duración de la Póliza", 
       x = "Duración de la Póliza (días)", 
       y = "Densidad") +
  theme(text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12),
        panel.border = element_blank(),
        plot.margin = margin(1, 1, 1, 1, "cm")) 

La distribución de la duración de las pólizas muestra un comportamiento relativamente uniforme a lo largo del eje x, con el máximo alrededor de los 175 días. Los valores mínimos se observan al inicio y hacia el final del rango, indicando una menor frecuencia de pólizas con duraciones extremadamente cortas o largas.

Mostrar código
set.seed(1)
uniforme <- runif(n = nrow(train), min = min(train$duracion_poliza), max = max(train$duracion_poliza))

data_comparacion <- data.frame(
  duracion_poliza = c(train$duracion_poliza, uniforme),
  tipo = rep(c("Real", "Uniforme"), each = nrow(train))
)
  
ggplot(data_comparacion, aes(x = duracion_poliza, fill = tipo, color = tipo)) +
  geom_density(alpha = 0.6) +
  theme_minimal() +
  labs(title = "Comparación de Distribución Duración Póliza y Uniforme", 
       x = "Duración de la Póliza (días)", 
       y = "Densidad") +
  scale_fill_manual(values = c(colores[1], colores2[1])) +
  scale_color_manual(values = c(colores[5], colores2[5])) +
  theme(text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12),
        legend.title = element_blank(),
        panel.border = element_blank(),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  guides(fill = guide_legend(override.aes = list(alpha = 0.6)))

Simulando una distribución uniforme aleatoria (con la semilla 1) podemos ver claramente que son en efecto bastante similares. Esto implica que todas las duraciones tienen una probabilidad similar de ocurrir, lo cual es consistente con una distribución uniforme, donde cada valor dentro del rango tiene la misma probabilidad.

Además, puede ser interesante estudiar el límite máximo de la variable:

max(train$duracion_poliza)
[1] 366

Que el valor máximo de la variable sea 366 días sugiere que el contrato abarca un máximo de un año calendario completo en un año bisiesto, aunque no es lo más común.

cor(train$duracion_poliza, train$numerosi)
[1] 0.2735884

Para interpretar los valores de correlación vamos a seguir estos intervalos:

Intervalos de interpretación del coeficiente de correlación:
  • \(|r|=1\): Correlación perfecta.
  • \(|r|\geqslant0.8\): Correlación alta.
  • \(|r|\geqslant0.5\): Correlación media
  • \(|r|\geqslant0.3\): Correlación débil
  • \(|r|>0\): Correlación muy débil
  • \(r=0\): No hay correlación.

La correlación entre la duración de la póliza y el número de siniestros es muy débil pero positiva (\(r=0.27\)). Esto sugiere que, en general, a medida que aumenta la duración de la póliza, también tiende a incrementarse el número de siniestros, aunque la relación no es especialmente fuerte. Una póliza más larga brinda más tiempo durante el cual puede ocurrir un siniestro. Por tanto, el aumento en el número de siniestros es razonable desde una perspectiva temporal.

De todas formas, aunque la relación existe, no es determinante. Otros factores además de la duración probablemente influyen significativamente en el número de siniestros, como el tipo de póliza, la antigüedad del vehículo o las características del asegurado. Podría ser interesante estudiar la media de siniestros para cada cantidad de días de duración.

Mostrar código
df_plot <- data.frame(
  duracion_poliza = sort(unique(train$duracion_poliza)),
  promedio_siniestros = tapply(train$numerosi, train$duracion_poliza, mean)
)

ggplot(df_plot, aes(x = duracion_poliza, y = promedio_siniestros)) +
  geom_line(color = colores[1], size = 1.2) +
  geom_point(color = colores[5], size = 3, shape = 19) + 
  theme_minimal() +
  labs(title = "Promedio de Siniestros por Duración de la Póliza",
       x = "Duración de la Póliza (días)", 
       y = "Promedio de Siniestros") +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"), 
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    axis.text = element_text(size = 10), 
    panel.grid.major = element_line(color = "gray", size = 0.5), 
    panel.grid.minor = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")  
  )

Cada punto en el gráfico muestra el promedio de siniestros para una duración específica de la póliza. Parece que hay más variabilidad en el promedio de siniestros a medida que aumenta la duración de la póliza. Esto significa que, conforme pasa el tiempo, es más común que las pólizas tengan diferentes números de reclamaciones. No se observa una hay una tendencia muy clara en los puntos del gráfico, lo que sugiere que el número promedio de siniestros no tiene un patrón fácil de predecir solo con la duración de la póliza.

3.2.2 Potencia del vehículo

summary(train$potencia)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  35.00   60.00   80.00   86.75  103.00  272.00 
Mostrar código
ggplot(train, aes(x = potencia)) +
  geom_density(fill = colores[1], color = colores[5], alpha = 0.6) + 
  theme_minimal() +
  labs(title = "Distribución de la Densidad de la Potencia", 
       x = "Potencia", 
       y = "Densidad") +
  theme(text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12),
        panel.border = element_blank(),
        plot.margin = margin(1, 1, 1, 1, "cm")) 

Esta variable parece seguir una distribución similar a la log-normal. Lo comprobamos:

Mostrar código
set.seed(1)
log_normal <- rlnorm(n = nrow(train), 
                     meanlog = mean(log(train$potencia + 1)), 
                     sdlog = sd(log(train$potencia + 1)))

data_comparacion_potencia <- data.frame(
  potencia = c(train$potencia, log_normal),
  tipo = rep(c("Real", "Log-Normal"), each = nrow(train))
)

ggplot() +
  geom_density(data = subset(data_comparacion_potencia, tipo == "Real"), 
               aes(x = potencia, fill = "Real"), 
               alpha = 0.6, color = colores[5]) +
  geom_density(data = subset(data_comparacion_potencia, tipo == "Log-Normal"), 
               aes(x = potencia, fill = "Log-Normal"), 
               alpha = 0.4, color = colores2[5]) +
  theme_minimal() +
  labs(title = "Comparación de Distribución Potencia y Log-Normal", 
       x = "Potencia", 
       y = "Densidad") +
  scale_fill_manual(values = c("Real" = colores[1], "Log-Normal" = colores2[1])) +
  theme(text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12),
        legend.title = element_blank(),
        panel.border = element_blank(),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  guides(fill = guide_legend(override.aes = list(alpha = 0.6)))

La distribución real de la potencia y la simulada con una log-normal son muy similares, lo que sugiere que la potencia de los vehículos en el conjunto de datos podría aproximarse razonablemente a este tipo de distribución. Esto puede darse porque en muchas situaciones reales, las variables relacionadas con características físicas o de rendimiento (como la potencia del motor) tienden a distribuirse de forma asimétrica, con una mayoría de valores en un rango medio y una cola hacia valores más altos que en este caso serían los coches de carreras ode lujo. Este comportamiento es característico de una distribución log-normal.

Ahora estudiamoscómo varía el número promedio de siniestros (variable respuesta) en función de la potencia del vehículo:

medias <- tapply(train$numerosi, train$potencia, mean)
Mostrar código
plot(sort(unique(train$potencia)), medias, pch = 20, col = colores[3], 
     xlab = "Potencia del vehículo", 
     ylab = "Número promedio de siniestros", 
     main = "Relación entre Potencia y Número Promedio de Siniestros", 
     cex = 1.2,   # Tamaño de los puntos
     cex.axis = 1, # Tamaño de los ejes
     cex.lab = 1.1, # Tamaño de las etiquetas
     cex.main = 1.3, # Tamaño del título
     lwd = 2,      # Grosor de la línea de los puntos
     lty = 1,      # Tipo de línea (sin línea, solo puntos)
     xlim = c(min(sort(unique(train$potencia))) - 10, max(sort(unique(train$potencia))) + 10))

Aparentemente, la mayoría de los puntos están agrupados en los valores más bajos de potencia y número promedio de siniestros, lo que sugiere que los vehículos con menor potencia tienden a tener menos siniestros. A parte, hay algunos puntos fuera de la masa principal de puntos, con valores más altos tanto en potencia como en el número promedio de siniestros. Esto indica que, aunque no es la norma, algunos vehículos con alta potencia tienen un número significativo de siniestros.

En nuestra base de datos no parece haber una tendencia específica, ya que los puntos se distribuyen al rededor de 0. Esto podría sugerir que los vehículos de menor potencia son más seguros o que son conducidos de manera más conservadora.

cor(train$potencia, train$numerosi)
[1] 0.08484551

Volviendo a la interpretación de cada rango de correlación, un valor de \(r=0.08\) indica una relación lineal directa muy débil, lo cual podíamos esperar del gráfico anterior ya que no se observa ninguna relación clara.

La discretización de una variable continua, como en este caso la potencia de un vehículo, implica dividirla en grupos o categorías. Este proceso puede ser útil para facilitar la interpretación de la relación entre esa variable y otras, como el número de siniestros. Se pueden usar los cuantiles de la variable para dividir los datos en grupos más equitativos en términos de cantidad de observaciones.

# Discretizar la potencia usando cuantiles
qnt <- quantile(train$potencia, seq(0, 1, .2))
qnt[1] <- qnt[1] - 0.01  # Ajustar el primer cuantíl
qnt[5] <- qnt[5] + 0.01  # Ajustar el quinto cuantíl

train$g.potencia <- cut(train$potencia, breaks = qnt)

train$g.potencia <- factor(train$g.potencia, levels = levels(train$g.potencia))

medias <- tapply(train$numerosi, train$g.potencia, mean)
Mostrar código
ggplot(data = data.frame(x = 1:length(medias), y = medias), aes(x = x, y = y)) +
  geom_point(color = colores[5], size = 2) + 
  geom_line(color = colores[3], size = 0.5) + 
  scale_x_continuous(breaks = 1:length(medias), labels = levels(train$g.potencia)) +
  labs(title = "Promedio de Siniestros por Grupo de Potencia", 
       x = "Rango de Potencia", 
       y = "Promedio de Siniestros") + 
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

Como podemos ver, los vehículos en el rango de potencia (90-112) parecen estar asociados con un mayor número de siniestros. Esto podría deberse a que estos vehículos tienen suficiente potencia para ser utilizados de manera más agresiva o en condiciones que aumentan el riesgo de siniestros.

Los vehículos en el rango de potencia (35-90) y (112-272) tienen menos siniestros en promedio, lo cual podría sugerir que estos vehículos son más seguros o que son conducidos de manera más conservadora, de hecho probablemente estos vehículos son tractores, patinetes y otros similares.

La variabilidad en otros rangos muestra que hay diferencias aparentemente significativas en el número promedio de siniestros según la potencia del vehículo.

La discretización de la variable ha sido muy interesante, podemos estar interesados en hacerlo de nuevo, por eso hemos creado una función para hacerlo de forma más automática:

Mostrar código
plot.dist <- function(predictor, cortes=8, variable=F){
      qnt <- quantile(predictor, seq(0,1,1/cortes))
      qnt[1] <- qnt[1]-0.01; qnt[cortes] <- qnt[cortes]+0.01
      V.discreta <- cut (predictor,  breaks= qnt)
      medias <- tapply(train$numerosi, V.discreta, mean)
      if (variable == T) return(V.discreta)
} # Fin función

Esta función nos ayudará a realizar los cálculos de discretización de forma mucho más rápida.

3.2.3 Antigüedad del vehículo

Estudiamos ahora la variable de antigüedad del vehículo.

Mostrar código
ggplot(train, aes(x = antivehi)) +
  geom_density(bw = 0.5, fill = colores[3], color = colores[5], alpha = 0.6) +
  labs(
    title = "Distribución de la Antigüedad del Vehículo",
    x = "Antigüedad del Vehículo (años)",
    y = "Densidad"
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

El gráfico no parece ajustarse a una distribución específica, pero sí muestra una clara tendencia decreciente. Observamos que hay muchos vehículos relativamente nuevos, con menos de 5 años de antigüedad, seguido de una disminución significativa a partir de entre los 11 y 12 años.

Para realizar la comparación, identificamos que la distribución más similar a nuestros datos es la chi-cuadrado, aunque no parece tan similar como en los casos anteriores. Esta distribución se caracteriza por su asimetría, comenzando con valores altos que disminuyen progresivamente. Es adecuada en este caso porque refleja una tendencia inicial elevada que va decreciendo, algo similar a lo observado en nuestros datos de antigüedad.

Mostrar código
scale_chi <- mean(train$antivehi) / 2

set.seed(1)
chi_data <- rchisq(n = length(train$antivehi), df = 2) * scale_chi

data_comparacion_antivehi <- data.frame(
  antiguedad = c(train$antivehi, chi_data),
  tipo = c(rep("Real", length(train$antivehi)), rep("Chi-Cuadrado", length(chi_data)))
)


ggplot() +
  geom_density(data = subset(data_comparacion_antivehi, tipo == "Real"), 
               aes(x = antiguedad, fill = "Real"), 
               alpha = 0.6, color = colores[5]) +
  geom_density(data = subset(data_comparacion_antivehi, tipo == "Chi-Cuadrado"), 
               aes(x = antiguedad, fill = "Chi-Cuadrado"), 
               alpha = 0.4, color = colores2[5]) +
  theme_minimal() +
  labs(title = "Comparación de Distribución Antigüedad del Vehículo y Chi-Cuadrado", 
       x = "Antigüedad del Vehículo (años)", 
       y = "Densidad") +
  scale_fill_manual(values = c("Real" = colores[1], "Chi-Cuadrado" = colores2[1])) +
  theme(text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12),
        legend.title = element_blank(),
        panel.border = element_blank(),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  guides(fill = guide_legend(override.aes = list(alpha = 0.6)))

Como podemos observar, ambas distribuciones son parecidas, salvo por un repunte en los datos reales alrededor de los 10 años de antigüedad. Este comportamiento podría estar relacionado con un mayor número de vehículos que alcanzan esa antigüedad debido a características del mercado o patrones de uso. A pesar de esta diferencia, la similitud se recupera al volver a descender posteriormente.

Mostrar código
antivehi_data <- data.frame(
  antivehi = sort(unique(train$antivehi)),
  mean_numerosi = tapply(train$numerosi, train$antivehi, mean)
)

ggplot(antivehi_data, aes(x = antivehi, y = mean_numerosi)) +
  geom_point(color = colores[5], size = 2) + 
  geom_line(color = colores[3], size = 0.5) + 
  theme_minimal() +
  labs(
    title = "Relación entre Antigüedad del Vehículo y Número de Siniestros",
    x = "Antigüedad del Vehículo (años)",
    y = "Promedio de Número de Siniestros"
  ) +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

Los vehículos más nuevos (0 años) tienen el promedio de siniestros más alto, más de 0.87. Esto podría deberse a la falta de experiencia del conductor con el vehículo nuevo. A medida que los vehículos son ligeramente más antiguos (de 1 a 5 años), el promedio de siniestros disminuye, alcanzando un mínimo local a los 7 años. Después, el promedio de siniestros muestra algunas fluctuaciones, pero en general vemos que decrece con algunos incrementos ocasionales.

Tras ver el gráfico, puede ser interesante ver la correlación entre la variable de respuesta, número de siniestros, y la antigüedad del vehículo.

cor(train$antivehi, train$numerosi) 
[1] -0.1807625

Hay una muy débil correlación directa (\(r=-0.18\)), es decir, existe una pequeña tendencia a que, a medida que aumenta la antigüedad del vehículo, disminuye el número de siniestros en promedio. El hecho de que sea tan débil implica que la antigüedad del vehículo no es un predictor particularmente fuerte del número de siniestros.

3.2.4 Tipo de vehículo

La variable tipo de vehículo es de tipo factor, por lo tanto es preciso otro tipo de análisis.

summary(train$tipovehi)
100 120 150 200 250 300 310 
581   8  32  31   2   3  10 

Estudiamos su estructura con un gráfico:

Mostrar código
ggplot(train, aes(x = tipovehi, y = numerosi)) +
  stat_summary(fun = mean, geom = "bar", fill = colores[1], color = colores[5], width = 0.6) +
  theme_minimal() +
  labs(title = "Promedio de Siniestros por Tipo de Vehículo", 
       x = "Tipo de Vehículo", 
       y = "Promedio de Siniestros") +
  theme(text = element_text(size = 12),
        plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        axis.title.x = element_text(size = 12), 
        axis.title.y = element_text(size = 12),
        axis.text.x = element_text(angle = 45, hjust = 1), 
        panel.border = element_blank(),
        plot.margin = margin(2, 1, 1, 1, "cm"))

El gráfico muestra el promedio de siniestros (reclamaciones) para diferentes tipos de vehículos. Los vehículos del tipo 120 y 310 tienen los promedios de siniestros más altos (aproximadamente 1 y 0.65 respectivamente). Esto sugiere que estos tipos de vehículos tienen más probabilidades de estar involucrados en siniestros. Los vehículos del tipo 250 y 300 no reportan a penas siniestros y los del resto de tipos tienen un promedio bajo (menos de 0.5). Esto indica que estos vehículos son menos propensos a involucrarse en siniestros.

Hay diferencias significativas en el promedio de siniestros entre los diferentes tipos de vehículos, lo cual es importante para los estudios de evaluación de riesgos en seguros. Estas diferencias hacen que la distribución no se asemeje a ninguna distribución común.

Ahora, al ser la variable de tipo factor, comprobamos si hay diferencias significativas en el número de siniestros entre los distintos tipos de vehículos. Usando ANOVA sabemos si tipo de vehículo tiene un impacto importante en el número de siniestros o si es más probable que las diferencias que vemos sean solo por casualidad.

mod1 <- lm(train$numerosi ~ train$tipovehi)

anova(mod1)
Analysis of Variance Table

Response: train$numerosi
                Df Sum Sq Mean Sq F value Pr(>F)
train$tipovehi   6   4.53 0.75521  0.8034 0.5675
Residuals      660 620.43 0.94005               
Umbrales comunes de significancia

El coeficiente de correlación mide la intensidad y dirección de la relación lineal entre dos variables:

  • \(p-value=0.05\): Este es el umbral más comúnmente utilizado.
  • \(p-value=0.01\): Un umbral más estricto.
  • \(p-value=0.1\): A veces se usa en investigaciones exploratorias, pero no se considera concluyente.

El análisis de varianza muestra que el tipo de vehículo no tiene un impacto significativo en el número de siniestros. El valor F es bajo y el valor p es alto (0.5675), no entra en ninguno de los casos mencionados, lo que sugiere que las diferencias en el número de siniestros entre los diferentes tipos de vehículos podrían ser atribuibles al azar y no son estadísticamente significativas. En otras palabras, el tipo de vehículo no parece ser un factor importante para predecir el número de siniestros en este caso.

confint(mod1)
                       2.5 %    97.5 %
(Intercept)        0.3754061 0.5333719
train$tipovehi120 -0.1321019 1.2233239
train$tipovehi150 -0.4563302 0.2350522
train$tipovehi200 -0.4827444 0.2191277
train$tipovehi250 -1.8028934 0.8941154
train$tipovehi300 -1.5563821 0.6476041
train$tipovehi310 -0.4615821 0.7528041

Aquí vemos los intervalos de confianza al 95% para los coeficientes estimados del modelo lineal que hemos creado. Específicamente, estos números nos dicen el rango dentro del cual creemos que se encuentra el valor real del coeficiente, con un 95% de certeza.

Para todos los tipos de vehículo mencionados, los intervalos de confianza incluyen el cero en su rango, lo que indica que, en efecto, no podemos afirmar que exista una relación significativa entre estos tipos de vehículos y el número de siniestros. En otras palabras, no podemos concluir que el tipo de vehículo tenga un impacto claro sobre el número de siniestros en el modelo.

3.2.5 Plazas de vehículo

Estudiamos ahora la variable plazas de vehículo.

table(train$plazasvehi)

  5   6 
662   5 

Como vemos, es una variable muy desbalanceada, aunque en términos generales, plazasvehi tiene dos categorías (5 y 6 plazas), la gran mayoría de los vehículos tiene 5 plazas, mientras que solo una pequeña proporción tiene 6 plazas.

tapply(train$numerosi, train$plazasvehi, mean)
        5         6 
0.4471299 0.6000000 

Esta funciónnos permite ver cómo varía el promedio de siniestros (numerosi) según el número de plazas del vehículo (plazasvehi), lo cual elimina el efecto del desbalance. El resultado nos indica que para los vehículos con 5 plazas, el promedio de los siniestros (train$numerosi) es 0.4471, y para los vehículos con 6 plazas, el promedio de los siniestros es 0.6.

Los vehículos con 6 plazas suelen ser más grandes, lo que podría implicar que son utilizados de manera diferente o para viajes más largos, lo que podría aumentar las probabilidades de que ocurran siniestros. En cambio, los vehículos más pequeños de 5 plazas a menudo son conducidos de manera diferente o menos veces, lo que podría reducir la exposición al riesgo de siniestros.

3.2.6 Edad

Ahora estudiamos la variable de edad. Realizamos de nuevo el gráfico de la distribución de la variable:

ggplot(train, aes(x = edad)) +
  geom_density(fill = colores[3], color = colores[5], alpha = 0.6, bw = 0.5) +
  theme_minimal() +
  labs(title = "Distribución de la Edad", 
       x = "Edad", 
       y = "Densidad") +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

La distribución de la edad sigue una forma que se asemeja bastante a una distribución normal, lo cual es común en los conjuntos de datos con esta variable. Las vemos constrastadas en el siguiente gráfico:

Mostrar código
media_edad <- mean(train$edad, na.rm = TRUE)
sd_edad <- sd(train$edad, na.rm = TRUE)

ggplot(train, aes(x = edad)) +
  geom_density(aes(fill = "Datos Reales"), color = colores[5], alpha = 0.6, bw = 0.5) + 
  stat_function(fun = dnorm, args = list(mean = media_edad, sd = sd_edad), 
                color = colores2[5], size = 0.6,
                geom = "area", aes(fill = "Normal Ajustada"), alpha = 0.3) + 
  theme_minimal() +
  labs(title = "Comparación de Distribución Normal y Edad", 
       x = "Edad", 
       y = "Densidad") +
  scale_fill_manual(name = "Distribuciones", 
                    values = c("Datos Reales" = colores[3], "Normal Ajustada" = colores2[3])) + 
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

La distribución observada en los datos reales sigue una forma similar a una distribución normal, pero presenta algunas irregularidades. El gráfico muestra un pico notable cerca de los 50 años, lo que indica que la mayoría de los casos se concentran en este rango de edad. Además, se observa un descenso en la densidad justo después de los 50 años, lo que sugiere una disminución en la frecuencia de los datos en ese grupo, para el cual no queda muy clara la interpretación.

Estudiamos ahora su relación con la variable respuesta:

cor(train$edad, train$numerosi) 
[1] -0.008924572

Para empezar, el coeficiente de correlación ($r=-0.008) nos indica una correlación muy débil negativa. A medida que la edad del conductor aumenta, el número de siniestros disminuye en sentido general, pero es tan débil la relación que no nos da ninguna información relevante.

Ahora vemos esta relación graficada:

Mostrar código
edad_data <- data.frame(
  edad = sort(unique(train$edad)),
  mean_numerosi = tapply(train$numerosi, train$edad, mean)
)

ggplot(edad_data, aes(x = edad, y = mean_numerosi)) +
  geom_point(color = colores[5], size = 2) + 
  geom_line(color = colores[3], size = 0.5) + 
  theme_minimal() +
  labs(
    title = "Relación entre Edad y Número de Siniestros",
    x = "Edad",
    y = "Promedio de Número de Siniestros"
  ) +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

En efecto, no queda muy clara la relación entre las variables. Por lo tanto, tal como hemos aprendido, discretizar la variable puede ayudar a su interpretación. La discretizamos usando la función que habíamos creado previamente.

dist.test <- function(predictor.train, predictor.test, cortes=10) {
    qnt <- quantile(predictor.train, seq(0, 1, 1/cortes))
    qnt[1] <- qnt[1] - 0.01
    qnt[cortes] <- qnt[cortes] + 0.01
    V.discreta <- cut(predictor.test, breaks = qnt)
    return(V.discreta)
}

train$edad_discretizada <- dist.test(train$edad, train$edad, cortes = 10)
test$edad_discretizada <- dist.test(train$edad, test$edad, cortes = 10)
levels(test$edad_discretizada) <- levels(train$edad_discretizada)
Mostrar código
edad_data_discretizada <- data.frame(
  edad_discretizada = levels(train$edad_discretizada),
  mean_numerosi = tapply(train$numerosi, train$edad_discretizada, mean)
)

edad_data_discretizada$edad_discretizada <- factor(edad_data_discretizada$edad_discretizada, 
                                                   levels = edad_data_discretizada$edad_discretizada)

ggplot(edad_data_discretizada, aes(x = as.numeric(edad_discretizada), y = mean_numerosi)) +
  geom_point(color = colores[5], size = 2) + 
  geom_line(color = colores[3], size = 0.5, aes(group = 1)) +  
  scale_x_continuous(breaks = 1:length(edad_data_discretizada$edad_discretizada), 
                     labels = edad_data_discretizada$edad_discretizada) + 
  theme_minimal() +
  labs(
    title = "Relación entre Edad Discretizada y Promedio de Número de Siniestros",
    x = "Grupo de Edad Discretizado",
    y = "Promedio de Número de Siniestros"
  ) +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

Ahora la relación queda mucho más clara. Los conductores más jóvenes (36-40 años) tienen un promedio de siniestros bastante alto, probablemente debido a la falta de experiencia al volante y comportamientos de conducción más arriesgados. Los conductores de edad joven-media (20-36) muestra una disminución en el promedio de número de siniestros, lo cual indica un mejor manejo al volante o hábitos más seguros al volante con los años. En los conductores de edad media (40-53) disminuye el promedio de siniestros, lo cual puede deberse a que los conductores en este grupo de edad, además de tener más experiencia al volante, a menudo tienen una vida profesional y familiar activa, por lo que podría significar que conducen con más cuidado. Los conductores más mayores (53-84) tienen un promedio de siniestros más bajo, posiblemente porque conducen con mayor precaución o debido a que no suelen conducir largos y complicados viajes sino cortos viajes del día a día para lo esencial.

3.2.7 Sexo conductor

En este caso analizamos la variable del sexo del conductor. Esta variable vuelve a ser de tipo factor así que la vamos a analizar de forma parecida a la variable de tipo ed vehículo.

table(train$sexocondu)

  M   V 
131 536 

De acuerdo con los datos, la mayoría de las conductoras son femeninas (dando por hecho que M es masculino y V femenino), ya que hay 536 observaciones de conductores femeninos en comparación con 131 de conductores masculinos. Esto indica que el grupo de conductores en este conjunto de datos está desbalanceado en favor de las mujeres.

tapply(train$numerosi, train$sexocondu, mean)
        M         V 
0.5648855 0.4197761 

Como vemos, los hombres en esta base de datos tienden a tener más siniestros en comparación con las mujeres, pero debemos comprobar si hay diferencias significativas en el número de siniestros entre estos. Usando ANOVA sabemos si el sexo tiene un impacto importante en el número de siniestros o si es más probable que las diferencias que vemos sean solo por casualidad.

mod2 <- lm(train$numerosi ~ train$sexocondu)

anova(mod2)
Analysis of Variance Table

Response: train$numerosi
                 Df Sum Sq Mean Sq F value Pr(>F)
train$sexocondu   1   2.22 2.21667  2.3671 0.1244
Residuals       665 622.75 0.93646               

Recordando los umbrales comunes de significancia, el análisis de varianza muestra que el sexo tiene un impacto significativo al nivel de significatividad 0.01 en el número de siniestros. El valor p es bajo (0.1244), lo que sugiere que las diferencias en el número de siniestros entre los dos sexos podrían deberse. La variable sexo del conductor parece ser un factor importante para predecir el número de siniestros.

confint(mod2)
                      2.5 %     97.5 %
(Intercept)       0.3988698 0.73090123
train$sexoconduV -0.3303046 0.04008584

Este intervalo de confianza sugiere que con un 95% de confianza, el efecto de ser mujer en la variable objetivo es negativo. Es decir, el ser mujer reduce el número de siniestros en comparación con los hombres tal como esperábamos con los estudios anteriores.

3.2.8 Antigüedad en la compañía

Ahora estudiamos la variable antigüedad en la compañía:

Mostrar código
ggplot(train, aes(x = anticia)) +
  geom_density(fill = colores[3], color = colores[5], alpha = 0.6, bw = 0.5) +
  theme_minimal() +
  labs(title = "Distribución de la Antigüedad en la Compañía",
       x = "Antigüedad del Vehículo (años)",
       y = "Densidad") +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

La mayoría de los clientes en la muestra llevan entre 0 y 5 años en la compañía Este rango tiene el pico de densidad más alto, cerca de 0.20, lo que sugiere que los clientes nuevos son los más comunes. A medida que la antigüedad aumenta, la densidad disminuye gradualmente. Esta distribución nos recuerda a una de tipo exponencial. Las vemos comparadas a continuación:

Mostrar código
lambda_anticia <- 1 / mean(train$anticia, na.rm = TRUE)

ggplot(train, aes(x = anticia)) +
  geom_density(aes(fill = "Datos Reales"), color = colores[5], alpha = 0.6, bw = 0.5) +  # Densidad observada
  stat_function(fun = dexp, args = list(rate = lambda_anticia), 
                color = colores2[5], size = 0.6,
                geom = "area", aes(fill = "Exponencial Ajustada"), alpha = 0.3) +  # Curva exponencial ajustada
  theme_minimal() +
  labs(title = "Comparación Distribución Exponencial y Antigüedad del Vehículo", 
       x = "Antigüedad del Vehículo (años)", 
       y = "Densidad") +
  scale_fill_manual(name = "Distribuciones", 
                    values = c("Datos Reales" = colores[3], "Exponencial Ajustada" = colores2[3])) + 
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

Como podemos observar, las distribuciones son bastante similares. Sin embargo, la distribución real no comienza tan alta como la exponencial y presenta un descenso más pronunciado entre aproximadamente los 3 y los 7 años. Además, hacia el final, la distribución real no muestra una caída tan pronunciada, sino más mantenida. Aunque las dos distribuciones son similares, las diferencias en el inicio, el comportamiento intermedio y el final indican que la distribución real tiene características particulares que no se ajustan perfectamente al modelo exponencial, sugiriendo que otros factores están influyendo en la antigüedad de los clientes en la compañía.

Ahora estudiamos su relación con la variable respuesta:

cor(train$anticia, train$numerosi)
[1] -0.01805373

Observamos una correlación muy débil (\(r=-0.01\)), lo que sugiere que, en general, a medida que un cliente lleva más años con la empresa, tiende a tener menos accidentes con el vehículo. Sin embargo, esta relación es tan débil que no es lo suficientemente “significativa” como para considerarse relevante.

Mostrar código
anticia_data <- data.frame(
  anticia = sort(unique(train$anticia)),
  mean_numerosi = tapply(train$numerosi, train$anticia, mean)
)

ggplot(anticia_data, aes(x = anticia, y = mean_numerosi)) +
  geom_point(color = colores[5], size = 2) +  # Puntos en el gráfico
  geom_line(color = colores[3], size = 0.7, aes(group = 1)) +  # Línea de conexión entre puntos
  labs(
    title = "Relación entre Antigüedad del Vehículo y Promedio de Siniestros",
    x = "Antigüedad del Vehículo (años)",
    y = "Promedio de Número de Siniestros"
  ) +
  theme_minimal() +  # Estilo minimalista
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    panel.border = element_blank(),
    plot.margin = margin(1, 1, 1, 1, "cm")
  )

De nuevo, el gráfico no deja clara ninguna relación entre las variables, así que discretizamos la variable agrupándola en si lleva menos de 3 años en la empresa (nuevo o new) o más (veterano o vet).

train$g.anticia <- ifelse(train$anticia < 3, "New", "Vet") 

train$g.anticia <- as.factor(train$g.anticia)
tapply(train$numerosi, train$g.anticia, mean)
      New       Vet 
0.5327635 0.3544304 

A pesar de que los vehículos nuevos tienen un promedio ligeramente mayor de siniestros que los vehículos veteranos, la diferencia es bastante pequeña. Esto sugiere que la antigüedad del vehículo, según esta discretización en “New” y “Vet”, no parece tener un impacto significativo en el número promedio de siniestros.

3.2.9 Usos del Vehículo

Por último, analizamos la variable uso del vehículo que ya habíamos discretizado en dos grupos al inicio de la práctica.

table(train$usovehi2)

 G1  G2 
627  40 
tapply(train$numerosi, train$usovehi2, mean)
       G1        G2 
0.4529506 0.3750000 

Los vehículos con un uso más alto (grupo G1 con uso mayor a 500) tienen un promedio de siniestros más alto en comparación con los vehículos con un uso menor (grupo G2 con uso menor a 100). Esto podría sugerir que los vehículos con más uso (G1) están involucrados en más accidentes o siniestros en promedio.

La convertimos a tipo factor para hacer un análisis más profundo de estas diferencias:

train$usovehi2 <- as.factor(train$usovehi2)

mod3 <- lm(train$numerosi ~ train$usovehi2)

anova(mod3)
Analysis of Variance Table

Response: train$numerosi
                Df Sum Sq Mean Sq F value Pr(>F)
train$usovehi2   1   0.23 0.22848  0.2432 0.6221
Residuals      665 624.74 0.93945               

Dado que el valor p es alto, 0.6221, no hay evidencia estadística suficiente para concluir que el uso del vehículo tiene un efecto significativo sobre el número de siniestros. En otras palabras, las diferencias entre los grupos G1 y G2 en términos de su promedio de siniestros podrían darse simplemente por el azar.

confint(mod3)
                      2.5 %    97.5 %
(Intercept)       0.3769454 0.5289558
train$usovehi2G2 -0.3883180 0.2324169

El intervalo de confianza que incluye el valor 0 indica que la diferencia en el efecto entre el grupo G1 y el grupo G2 no es estadísticamente significativa. Por lo tanto, como ya hemos visto, no podemos afirmar que ser parte del grupo G2 tenga un impacto significativo en el número de siniestros.

3.3 Algunas conclusiones del EDA

En el análisis exploratorio de datos (EDA), la mayoría de las variables estudiadas y sus respectivas transformaciones parecen tener una capacidad predictiva notable por sí solas, con algunas excepciones como las variables plazas del vehículo y tipo de vehículo, que no muestran un impacto significativo. Por tanto, se optará por incluir todas estas variables en el modelo inicial.

En el caso de edad o antigüedad del carnet, hemos observado una relación no lineal con la variable dependiente. Este comportamiento sugiere que la edad podría no seguir una tendencia lineal simple, y es por eso que hemos decidido incluir transformaciones de la edad, como edad al cuadrado y edad al cubo, que pueden captar esta no linealidad de manera más efectiva, que funciona como alternativa al uso de los modelos GAM.

# Edad cuadrado y cubo
train$edad2 <- train$edad^2
train$edad3 <- train$edad^3

Hay muchas oportunidades para construir más predictores que podrían mejorar el modelo. Algunas posibles ideas incluyen:

  • Discretizar variables adicionales o hacerlo de manera diferente, por ejemplo, agrupando rangos de edad, antigüedad del vehículo, etc.
  • Incorporar efectos de interacción entre variables, ya que algunas combinaciones de variables podrían tener un impacto mayor al combinarse.
  • Crear nuevas variables derivadas a partir de las que ya tenemos disponibles. Un ejemplo sería utilizar el código postal para generar una nueva variable que capture la localización geográfica y su influencia potencial en el comportamiento de los siniestros o cualquier otra variable de interés.

Esto puede enriquecer el modelo y mejorar su capacidad predictiva, por lo que seguir explorando nuevas combinaciones y transformaciones de variables sería una parte clave del proceso de modelado en otro caso.

4 Evaluando modelos

Propuestas de modelización

Especificación:

  • Variables que van a entrar en el predictor lineal
  • Distribución para la variable respuesta
  • Función link
  • Y con enfoque ML, tipo de regularización/penalización

Distribuciones para la variable respuesta

Como modelos para la variable respuesta podemos testar:

  • Poisson
  • Binomial Negativa
  • ZI (inflada en cero) + Poisson
  • ZI (inflada en cero) + BN
  • Normal

Funciones link

Como funciones link:

  • Logaritmo
  • Identidad

Regularización

Penalización: * Net Elastic (Lasso, Ridge)

Modelos sin y con regularización
  • Modelos sin regularización: Estos modelos ajustan los datos de forma directa, maximizando el ajuste sin restricciones. Aunque pueden ofrecer un buen desempeño inicial, son más propensos al sobreajuste, capturando ruido en lugar de patrones generales.

  • Modelos con regularización: Incorporan penalizaciones que limitan la complejidad del modelo, favoreciendo soluciones más simples y generalizables. Ayudan a evitar el sobreajuste y suelen ser más robustos en escenarios reales..

4.1 Modelos sin regularización

Cuando utilizamos modelos sin regularización hemos de estar atentos a los problemas de multicolinealidad.

Antes de comenzar a modelizar hemos de eliminar algunas variables por que es lo razonable, lógico, desde el punto de vista de modelización y de computación.

Eliminamos codigopostal y ccaa ya que no aportan información relevante al modelo y no necesariamente tienen una relación significativa o directa con el número de siniestros, ya que ambas variable son identificadores geográficos y podrían afectar a la complejidad del modelo en lugar de mejorarlo. También eliminamos las variables iniciop y finalp, debido a que ya hemos calculado con anterioridad la diferencia entre estas para determinar la duración de la póliza, por lo que mantener estas variables sería redundante y podría dificultar la interpretación del modelo.

train <- subset(train, select = -c(codigopostal, ccaa,
                                   iniciop, finalp))

Ajuste modelos (clásicos) de Poisson

Po.log <- glm(numerosi ~ . , data = train, 
              family = poisson(link = "log"))

# Obtenemos los valores iniciales y nos aseguramos que no haya coeficientes negativos o NA
starting.values <- coef(Po.log)
starting.values[starting.values < 0] <- 0.0001
starting.values[is.na(starting.values)] <- 0.0001

# Ajuste modelo Poisson con enlace identidad usando los valores iniciales
Po.id <- glm(numerosi ~ . , data = train, start = starting.values,
             family = poisson(link = "identity"))

Ajuste modelo (clásico) Binomial Negativo

BN.log <- glm.nb(numerosi ~ . , data = train, link = log)
# starting.values <- coef(BN.log)
# starting.values[starting.values < 0] <- 0.0001
# BN.id <- glm.nb(numerosi ~ . , data = train,
#               start = starting.values, link = identity)

Ajustes de modelos (clásicos) inflados de ceros

# Modelo inflado de ceros con Poisson
ZI.Po <- zeroinfl(numerosi ~ ., data = train, dist = "poisson")

# Modelo inflado de ceros con Binomial Negativa
ZI.BN <- zeroinfl(numerosi ~ ., data = train, dist = "negbin")

Normal

La variable de respuesta, numerosi, representa el número de siniestros para cada observación, por lo que se trata de una variable discreta que toma valores enteros no negativos entre los cuales es común observar una gran cantidad de ceros.

Debido a esto, realizar un ajuste con modelo Normal no tiene sentido ya que la Normal es una distribución que admite valores negativos y podría generar errores significativos en la estimación de los parámetros y predicciones poco precisas. Por lo tanto, descartamos este modelo.

4.1.1 El rol del tiempo de exposición: Modelos con Offsets (y pesos)

A continuación, realizamos diversos modelos para analizar el número de siniestros, incluyendo el manejo del tiempo de exposición utilizando offsets y pesos.

Mostrar código
# Poisson con offset
Po.log.off <- glm(numerosi ~ . ,
                  data = subset(train, select = -c(duracion_poliza)),
                  offset = log(train$duracion_poliza / 365.25), # Convertimos la duración a años  
                  family = poisson(link = "log"))

# Poisson con pesos (corrigiendo `weight` a `weights`)
Po.log.wgt <- glm(numerosi ~ . ,
                  data = subset(train, select = -c(duracion_poliza)),
                  weights = train$duracion_poliza / 365.25,  
                  family = poisson(link = "log"))

El término offset nos permite incorporar la duración de la póliza como una variable que ajusta la tasa de siniestros por unidad de tiempo, garantizando que la duración de la póliza se considere de manera proporcional al análisis. En este caso, la duración se convierte a años, y su efecto se incluye de manera multiplicativa en el modelo de Poisson, ajustando las tasas de siniestros en función del tiempo de exposición al riesgo. Por otro lado, el modelo con pesos utiliza la duración de la póliza como un peso, asignando mayor relevancia a las observaciones con una mayor duración, y corrigiendo así el modelo para que tenga en cuenta de manera adecuada la variable de tiempo sin necesidad de usar un término offset explícitamente.

# Comparación de AIC's
AIC(Po.log.off)
[1] 1156.978
AIC(Po.log.wgt) * nrow(train) / sum(train$duracion_poliza / 365.25)
[1] 1580.097

Observamos que el modelo con offset presenta un AIC más bajo que el modelo con pesos. Esto sugiere que el primer modelo se ajusta mejor a los datos, además de ser más interpretable en términos de tasas de siniestros ajustadas por la duración de la póliza.

Mostrar código
# Poisson con offset (especificación alternativa)
of_var <- log(train$duracion_poliza / 365.25)
formula <- numerosi ~ potencia + antivehi +
           tipovehi + plazasvehi + edad  +
           sexocondu + anticia + prov + duracion_poliza +
           usovehi2 + g.potencia + offset(of_var) +
           g.anticia + edad2 + edad3 + edad_discretizada

Po.log.off2 <- glm(formula, data = train, family = poisson)

# Binomial Negativa con offset
BN.log.off <- glm.nb(formula, data = train, link = log)

# Modelos inflados de ceros, con offset
formula2 <- formula(numerosi ~ potencia +
             antivehi + tipovehi + plazasvehi + edad + 
             sexocondu + anticia + 
             prov + usovehi2 + g.potencia +
             offset(log(duracion_poliza/365.25)) +
             g.anticia + edad2 + edad3
           | potencia + antivehi + tipovehi +
             plazasvehi + edad + sexocondu +
             anticia + prov + usovehi2 + g.potencia +
             g.anticia + edad2 + edad3 + edad_discretizada)

ZI.Po.off <- zeroinfl(formula2, data = train, dist = "poisson")
ZI.BN.off <- zeroinfl(formula2, data = train, dist = "negbin")

Realizamos otros modelos utilizando una fórmula más detallada, especificando los predictores. La variable respuesta aparece como un offset explícito, lo cual aclara su contribución al modelo.

La Binomial Negativa permite manejar la sobredispersión que puede llegar a ocurrir con Poisson, y los modelos inflados de cero se utilizan cuando hay una gran cantidad de ceros en la variable dependiente, en este caso el número de siniestros, condición que no puede ser capturada adecuadamente por Poisson o Binomial Negativa.

4.1.2 Evaluación clásica de los modelos

Procedemos a comparar los diferentes modelos predictivos utilizando el criterio AIC.

AICs <- c(AIC(Po.log), AIC(Po.id), AIC(BN.log), AIC(ZI.Po), AIC(ZI.BN),
          AIC(Po.log.off), 
          # ‘Ajuste’ por número de observaciones
          AIC(Po.log.wgt)*nrow(train)/sum(train$duracion_poliza/365.25),
          AIC(BN.log.off), AIC(ZI.Po.off), AIC(ZI.BN.off))
Mostrar código
errores_clasicos <- data.frame(
  modelo = c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", 
              "Po.log.off", "Po.log.wgt", 
             "BN.log.off", "ZI.Po.off", "ZI.BN.off"), AICs)

errores_clasicos %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()
modelo AICs
Po.log 1178.052
Po.id 1279.821
BN.log 1149.058
ZI.Po 1141.794
ZI.BN 1117.741
Po.log.off 1156.978
Po.log.wgt 1580.097
BN.log.off 1131.617
ZI.Po.off 1098.081
ZI.BN.off 1126.908
Mostrar código
colores1 <- carto_pal(11, "Teal")

barplot(
  errores_clasicos$AICs,            
  names.arg = errores_clasicos$modelo,  
  col = colores1[rank(errores_clasicos$AICs)],    
  border = "white",                 
  main = "Valores AIC por Modelo",    
  xlab = "Modelos",                
  ylab = "AIC", 
  las = 1,                      
  cex.names = 0.6,
  cex.axis = 0.8, 
  cex.main = 1.2                 
)

Vemos que el modelo Poisson inflado de ceros con offset presenta el AIC más bajo (AIC = 1098.081), lo cual nos indica que es el modelo que mejor se ajusta a los datos. Elegiríamos este modelo, ya que un AIC bajo sugiere que este tiene una mejor relación entre la capacidad de ajuste y el número de parámetros, evitando así el sobreajuste.

4.1.3 Evaluación enfoque ML

Definición de nuevos predictores definidos en el conjunto de training en test.

Mostrar código
qnt <- quantile(train$potencia,seq(0,1,.2))

qnt[1] <- qnt[1]-0.01
qnt[5] <- qnt[5]+0.01

test$g.potencia <- cut(test$potencia,  breaks= qnt)

test$g.anticia <- ifelse(test$anticia < 3, "New", "Vet")

test$edad2 <- test$edad^2

test$edad3 <- test$edad^3

test <- subset(test, select = -c(codigopostal, ccaa, iniciop, finalp))
test <- test[test$duracion_poliza != 0,]

Discretizamos la potencia usando cuantiles en el test, como ya hicimos anteriormente en el training set. También clasificamos la antigüedad de la compañía en “Nueva” o “Veterana”, creamos tranformaciones para capturar relaciones no lineales entre edad y otras variables en modelos predictivos, eliminamos columnas innecesarias y por último, filtramos los datos para eliminar observaciones donde la duración de la póliza sea 0 para evitar errores posteriores.

Eliminación de observaciones

Tras haber intentado generar las predicciones, obtuvimos un error sobre los niveles 5 y 19 de la variable prov. En esta tabla observamos que dichos niveles no están presentes en el conjunto de entrenamiento mientras que sí lo están en el conjunto de prueba, generando un error ya que los modelos no sabrán tratar estos niveles. Para evitar errores durante la generación de predicciones y demás, eliminamos las observaciones de dichos niveles en el test set.

Mostrar código
resumen <- data.frame(
  Prov = names(table(train$prov)),
  Train = as.integer(table(train$prov)),
  Test = as.integer(table(test$prov)[names(table(train$prov))]),
  stringsAsFactors = FALSE
)

resumen[resumen$Train == 0,]
   Prov Train Test
5     5     0    1
19   19     0    1

Posteriormente, eliminamos del conjunto test los casos, para los que no se puede crear un predictor

test <- test[!(test$prov %in% c(5, 19)), ]
test0 <- test[rowSums(is.na(test))==0,]

Predicciones en el conjunto de test (modelos sin regularización)

Mostrar código
# Predicciones
p.Po.log <- predict(Po.log, test0 , type="response")
p.Po.id <- predict(Po.id, test0 , type="response")
p.BN.log <- predict(BN.log, test0 , type="response")
p.ZI.Po <- predict(ZI.Po, test0 , type="response")
p.ZI.BN <- predict(ZI.BN, test0 , type="response")
nd <- data.frame(test0, of_var = log(test0$duracion_poliza/365.25)) 
p.Po.log.off <- predict(Po.log.off2, nd , type = "response")
p.BN.log.off <- predict(BN.log.off, nd, type = "response")
p.ZI.Po.off <- predict(ZI.Po.off, test0 , type="response")
p.ZI.BN.off <- predict(ZI.BN.off, test0 , type="response")

Pasamos a generar las predicciones en el conjunto de prueba para los modelos sin regularización. Descartamos el modelo Poisson con pesos debido a que, como hemos podido observar, tiene un AIC significativamente más alto respecto al resto de modelos.

Mostrar código
colores3 <- carto_pal(11, "RedOr")

resultados <- data.frame(
  reales = test0$numerosi,
  Po_log = p.Po.log,
  Po_id = p.Po.id,
  BN_log = p.BN.log,
  ZI_Po = p.ZI.Po,
  ZI_BN = p.ZI.BN,
  Po_log_off = p.Po.log.off,
  BN_log_off = p.BN.log.off,
  ZI_Po_off = p.ZI.Po.off,
  ZI_BN_off = p.ZI.BN.off
)

resultados <- melt(resultados, id.vars = "reales", 
                   variable.name = "Modelo", 
                   value.name = "Predicciones")

ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +
  geom_density(alpha = 0.4) +
  geom_density(data = resultados, aes(x = reales), 
               fill = "black", alpha = 0.2) +
  facet_wrap(~ Modelo) +
  scale_x_continuous(limits = c(0, 5)) +
  scale_fill_manual(values = colores3) +
  labs(x = "Número de Siniestros", 
       y = "Densidad", 
       title = "Comparación de Predicciones con Datos Reales") +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    axis.text.x = element_text(),
    panel.border = element_blank(),
    plot.margin = margin(2, 1, 1, 1, "cm"),
    legend.position = "none")

El gráfico muestra la comparación de las distribuciones de las predicciones obtenidas con distintos modelos en el test set, frente a los datos reales del mismo. Observamos que en la mayoría de modelos las predicciones siguen una tendencia parecida a la de los datos reales, con un pico inicial elevado cerca de 0 sinestros. Mientras que modelos como Po_id presentan diferencias marcadas respecto a realidad, otros como ZI_BN_off parecen ajustarse mejor.

No obstante, para determinar correctamente cuál es el mejor modelo, procedemos a calcular los errores cuadráticos y absolutos.

Errores cuadráticos

Mostrar código
# Función de perdida cuadrática
e.Po.log <- sum((p.Po.log - test0$numerosi)^2)
e.Po.id <- sum((p.Po.id - test0$numerosi)^2)
e.BN.log <- sum((p.BN.log - test0$numerosi)^2)
e.ZI.Po <- sum((p.ZI.Po - test0$numerosi)^2)
e.ZI.BN <- sum((p.ZI.BN - test0$numerosi)^2)
e.Po.log.off <- sum((p.Po.log.off - test0$numerosi)^2)
e.BN.log.off <- sum((p.BN.log.off - test0$numerosi)^2)
e.ZI.Po.off <- sum((p.ZI.Po.off - test0$numerosi)^2)
e.ZI.BN.off <- sum((p.ZI.BN.off - test0$numerosi)^2)

e.base <- sum((mean(train$numerosi) - test0$numerosi)^2)
errores_cuadraticos <- c(e.Po.log, e.Po.id, e.BN.log, e.ZI.Po, e.ZI.BN,
                          e.Po.log.off, e.BN.log.off,
                         e.ZI.Po.off, e.ZI.BN.off, e.base)

El error cuadrático se utiliza para evaluar el ajuste de los modelos comparando las predicciones con los valores reales observados en el test set. Penaliza más fuertemente las grandes desviaciones entre las predicciones y los valores reales.

Errores valor absoluto

Mostrar código
a.Po.log <- sum(abs(p.Po.log - test0$numerosi))
a.Po.id <- sum(abs(p.Po.id - test0$numerosi))
a.BN.log <- sum(abs(p.BN.log - test0$numerosi))
a.ZI.Po <- sum(abs(p.ZI.Po - test0$numerosi))
a.ZI.BN <- sum(abs(p.ZI.BN - test0$numerosi))
a.Po.log.off <- sum(abs(p.Po.log.off - test0$numerosi))
a.BN.log.off <- sum(abs(p.BN.log.off - test0$numerosi))
a.ZI.Po.off <- sum(abs(p.ZI.Po.off - test0$numerosi))
a.ZI.BN.off <- sum(abs(p.ZI.BN.off - test0$numerosi))

a.base <- sum(abs(mean(train$numerosi) - test0$numerosi))
errores_abs <- c(a.Po.log, a.Po.id, a.BN.log, a.ZI.Po, a.ZI.BN, 
                a.Po.log.off, a.BN.log.off, 
                 a.ZI.Po.off, a.ZI.BN.off, a.base)

A diferencia del error cuadrático, el error absoluto no penaliza en exceso las grandes desviaciones, sino que simplemente calcula la suma de las diferencias absolutas entre las predicciones y los valores reales. Esto es útil cuando se quiere tener una idea más directa de la magnitud de los errores sin que los grandes errores afecten de forma desproporcionada el cálculo total.

Resumen errores (modelos sin regularización)

A continuación, se muestra una tabla resumiendo cuadráticos y absolutos obtenidos para cada modelo para facilitar la comparación entre los modelos ajustados.

Mostrar código
errores_sin <- data.frame(
  modelo = c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", 
              "Po.log.off", "BN.log.off", 
             "ZI.Po.off", "ZI.BN.off", "Null model"),
  errores_cuadraticos,
  errores_abs)

errores_sin[, 2:3] <- round(errores_sin[, 2:3], 2)
errores_sin %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()
modelo errores_cuadraticos errores_abs
Po.log 245.12 170.47
Po.id 257.66 199.57
BN.log 245.98 171.27
ZI.Po 914.22 209.36
ZI.BN 550.66 198.58
Po.log.off 268.97 174.67
BN.log.off 266.53 174.09
ZI.Po.off 382.59 192.51
ZI.BN.off 389.49 194.57
Null model 255.06 198.48
Mostrar código
errores_combinados <- rbind(
  errores_sin$errores_cuadraticos, 
  errores_sin$errores_abs
)

colnames(errores_combinados) <- c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", 
                                  "Po.log.off", "BN.log.off", "ZI.Po.off", "ZI.BN.off", "Null model")

barplot(errores_combinados, 
        beside = TRUE, 
        names.arg = c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", 
                      "Po.log.off", "BN.log.off", "ZI.Po.off", "ZI.BN.off", "Null model"), 
        col = c(colores2[2], colores2[7]),
        border = "white",                
        main = "Comparación de Errores Cuadráticos y Absolutos", 
        xlab = "Modelo", 
        ylab = "Error", 
        las = 1,                      
        cex.names = 0.6,
        cex.axis = 0.8, 
        cex.main = 1.2)

Atendiendo a los errores y a las características de los modelos, elegiríamos el modelo de la Binomial Negativa. Como se ha comentado anteriormente, la distribución Binomial Negativa se ajusta mejor a nuestra muestra y es útil a la hora de sobrellevar la sobredispersión en los datos que se genera comúnmente en variables donde la varianza supera la media. A pesar de que los errores más bajos los presenta el modelo Poisson con enlace log, los valores tanto del error cuadrático como del absoluto son muy cercanos entre ambos modelos.

Respecto a la evaluación clásica, aunque evaluar el AIC puede resultar en modelos más complejos con menos sobreajuste, los errores nos ofrecen una visión directa del rendimiento predictivo del modelo. Esto puede hacer que un modelo con un AIC más bajo pero con peores errores no será necesariamente la mejor opción si nos fijamos en la precisión de las predicciones.

4.2 Modelos con regularización

Tras una comparación inicial de modelos sin regularización basándonos en el ajuste de los datos, ahora incluiremos la regularización, penalizando la complejidad del modelo y promoviendo soluciones más simples y robustas.

Los Paquetes glmnet y mpath

Cuestiones sobre las funciones - Conversión de data.frame en matriz (construyendo variables dummy para factores e intercepto), necesario para predecir.

model.matrix(~., data)

  • Modelo con respuesta Poisson

cv.glmnet(x, y, family=“poisson”, alpha, nfolds) glmnet(x, y, family=“poisson”, alpha) cvglmreg(formula, data, family=“poisson”, penalty, nfolds) glmreg(formula, data, family=“poisson”, penalty)

  • Modelo con respuesta Binomial Negativa

cvglmregNB(formula, data, penalty, nfolds) glmregNB(formula, data, penalty)

  • Modelos inflados de ceros

cvzipath(formula, data, family, link, penalty, nfolds) zipath(formula, data, family, link, penalty)

  • Predicción

predict(object, newdata, type) predict(object, newx, which, type)

Más información

https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#poi https://cran.r-project.org/web/packages/mpath/mpath.pdf

Conversión a objetos matriz de predictores

# Convertimos los data.frame en objetos matriz respetando factores
train.m <- model.matrix(~., data=train)
train.m <- train.m[,colnames(train.m) != "numerosi"]

test.m <- model.matrix(~., data=test0)
test.m <- test.m[,colnames(test.m) != "numerosi"]

Modelos con regularización

Se estima un modelo, existe una solución, por cada valor del parámetro de penalización.

Descarte de modelos

A lo largo del resto del trabajo, hemos decidido descartar algunos modelos debido a errores técnicos que no se pudieron solucionar durante el proceso.

4.2.1 Ajustes con glmnet

Mostrar código
# Respuesta Normal
No.re <- glmnet(x = train.m, y = train$numerosi, alpha=0.5)

# Respuesta Poisson
Po.re <- glmnet(x = train.m, y = train$numerosi,
                family = "poisson", alpha = 0.5)

# Respuesta Poisson con offset
Po.re.off <- glmnet(x = train.m, y = train$numerosi, alpha = 0.5,  
                   offset = log(train$duracion_poliza/365.25), family = "poisson")

Utilizamos glmnet para ajustar modelos con regularización mediante la penalización Elastic Net (alpha = 0.5), que combina las propiedades de las penalizaciones Lasso y Ridge.

4.2.2 Ajustes con mpath

Mostrar código
# Respuesta Poisson
Po.re2 <- glmreg(numerosi ~ . , data = train,
                 family = "poisson", penalty="enet")

# Respuesta Binomial Negativa
BN.re <- glmregNB(numerosi ~ . , data = train, penalty="enet")

# Respuestas Infladas de Ceros
ZI.Po.re <- zipath(numerosi ~ .| . , data= train,
              family = "poisson", link = "logit", penalty="enet")

ZI.BN.re <- zipath(numerosi ~ .| . , data= train,
                   family = "negbin", link = "logit", penalty="enet")


# Modelos Po y BN con offset
formula3 <- numerosi ~ potencia + antivehi +
            tipovehi + plazasvehi + edad  +
            sexocondu + anticia + prov + usovehi2 +
            g.potencia + g.anticia + edad2 + edad3

# Modelos con offset
Po.re2.off <- glmreg(formula3, data = train, family = "poisson",
                 offset = log(duracion_poliza/365.25), penalty="enet")

BN.re.off <- glmregNB(formula3, data = train, 
                     offset = log(duracion_poliza/365.25), penalty="enet")

# Modelos ZI con offset
formula4 <- formula(numerosi ~ potencia +
                      antivehi + tipovehi + plazasvehi + edad  + sexocondu + anticia + 
                      prov + usovehi2 + g.potencia  + g.anticia + edad2 + edad3 | 
                      potencia + antivehi + tipovehi + plazasvehi + edad  + 
                      sexocondu + anticia + prov + usovehi2 + g.potencia  + g.anticia + edad2 + edad3)

ZI.Po.re.off <- zipath(formula4, data= train, family = "poisson", offset = log(duracion_poliza/365.25), link = "logit", penalty="enet", maxit = 100)

ZI.BN.re.off <- zipath(formula4, data= train, family = "negbin", offset = log(duracion_poliza/365.25), link = "logit", penalty="enet", maxit = 100)

Ajustamos los modelos con mpath con fórmulas específicas para asegurarnos de que las variables predictoras relevantes estén correctamente incluidas en los ajustes. Este enfoque nos permite comparar diferentes técnicas y distribuciones, considerando características particulares de los datos y explorando la influencia de la regularización en cada caso.

4.2.3 Elección del modelo

Disponemos de un modelo para cada lambda testado, ¿cuál elegimos?

Vamos a considerar dos alternativas:

  • Predecimos con el modelo con el lambda que genera el “mejor” ajuste con medida clásica (AIC).
  • Elegimos el mejor modelo por Cross-Validation

Elección de modelos con regularización basándonos en el AIC y predicción

Mostrar código
# Modelo de respuesta Normal
selec <- which(No.re$dev.ratio == max(No.re$dev.ratio))
p.No.re <- predict(No.re, test.m, type = "response")[, selec]

# Modelo de respuesta Poisson (glmnet)
selec <- which.max(Po.re$dev.ratio)
p.Po.re <- predict(Po.re, test.m, type = "response")[, selec]

# Modelo de respuesta Poisson y offset glmnet
selec <- which.max(Po.re.off$dev.ratio)
p.Po.re.off <- predict(Po.re.off, test.m, 
                       newoffset = log(test$duracion_poliza / 365.25),
                       type = "response")[, selec]

# Modelo de respuesta Poisson (mpath)
selec <- which.min(Po.re2$resdev)
p.Po.re2 <- predict(Po.re2, test0, which = selec, type = "response")

# Modelo de respuesta Poisson y offset (mpath)
selec <- which.min(Po.re2.off$resdev)
p.Po.re2.off <- exp(predict(Po.re2.off, test0, which = selec,
                            newoffset = log(test0$duracion_poliza / 365.25)))

# Modelo de respuesta BN (mpath)
selec <- which.min(AIC(BN.re))
p.BN.re <- predict(BN.re, test0, which = selec, type = "response")

# Modelo de respuesta BN y offset (mpath)
selec <- which.min(AIC(BN.re.off))
p.BN.re.off <- exp(predict(BN.re.off, test0, which = selec,
                           newoffset = log(test0$duracion_poliza / 365.25)))

# Modelo de respuesta ZI-Po (mpath)
selec <- which.min(AIC(ZI.Po.re))
p.ZI.Po.re <- predict(ZI.Po.re, test0, which = selec,
                      type = "response")

# Modelo de respuesta ZI-Po y offset (mpath)
selec <- which.min(AIC(ZI.Po.re.off))
p.ZI.Po.re.off <- predict(ZI.Po.re.off, test0, which = selec,
                          type = "response")

# Modelo de respuesta ZI-BN (mpath)
selec <- which.min(AIC(ZI.BN.re))
p.ZI.BN.re <- predict(ZI.BN.re, test0, which = selec, 
                      type = "response")

# Modelo de respuesta ZI-BN y offset (mpath)
selec <- which.min(AIC(ZI.BN.re.off))
p.ZI.BN.re.off <- predict(ZI.BN.re.off, test0, which = selec, 
                          type = "response")

Realizamos la selección y predicción de varios modelos ajustados con regularización, considerando qué tanto del comportamiento de los datos explica el modelo, o qué tan bueno es el modelo según el criterio del AIC. Además, incorporamos ajustes con offset para manejar correctamente la exposición al riesgo, asegurando que los modelos sean más interpretables y relevantes para el contexto.

Mostrar código
resultados <- data.frame(
  reales = test0$numerosi,
  No_re = p.No.re,
  Po_re = p.Po.re,
  Po_re_off = p.Po.re.off,
  Po_re2 = p.Po.re2,
  BN_re = p.BN.re,
  ZI_Po_re = p.ZI.Po.re,
  Po_re2_off = p.Po.re2.off,
  BN_re_off = p.BN.re.off,
  ZI_Po_re_off = p.ZI.Po.re.off,
  ZI_BN_re = p.ZI.BN.re,
  ZI_BN_re_off = p.ZI.BN.re.off
)

resultados <- melt(resultados, id.vars = "reales", 
                   variable.name = "Modelo", 
                   value.name = "Predicciones")

ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +
  geom_density(alpha = 0.4) +
  geom_density(data = resultados, aes(x = reales), 
               fill = "black", alpha = 0.2) +
  facet_wrap(~ Modelo) +
  scale_x_continuous(limits = c(0, 5)) +
  scale_fill_manual(values = colores3) +
  labs(x = "Número de Siniestros", 
       y = "Densidad", 
       title = "Comparación de Predicciones con Datos Reales") +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    axis.text.x = element_text(),
    panel.border = element_blank(),
    plot.margin = margin(2, 1, 1, 1, "cm"),
    legend.position = "none")

Comparamos las densidades de los predicciones de varios modelos con las de los datos reales del número de siniestros. Los modelos BN_re y BN_re_off destacan por ajustarse mejor a los datos reales, especialmente en los valores más frecuentes (0 y 1 siniestros). En cambio, modelos como Po_re_off tienden a sobreestimar los siniestros altos, lo que indica un ajuste menos preciso. Por otro lado, el modelo No_re, aunque tiene bajos errores, subestima frecuencias importantes y no representa bien la variabilidad de los datos. En general, los modelos basados en la Binomial Negativa logran capturar mejor las características reales de la distribución.

Ahora, al igual que hicimos antes con los modelos sin regularización, comprobamos los errores cuadráticos y absolutos.

Errores cuadráticos

Mostrar código
# Funcion de perdida cuadrática
e.No.re <- sum((p.No.re - test0$numerosi)^2)
e.Po.re <- sum((p.Po.re - test0$numerosi)^2)
e.Po.re.off <- sum((p.Po.re.off - test0$numerosi)^2)
e.Po.re2 <- sum((p.Po.re2 - test0$numerosi)^2)
e.BN.re <- sum((p.BN.re - test0$numerosi)^2)
e.ZI.Po.re <- sum((p.ZI.Po.re - test0$numerosi)^2)
e.ZI.BN.re <- sum((p.ZI.BN.re - test0$numerosi)^2)
e.Po.re2.off <- sum((p.Po.re2.off - test0$numerosi)^2)
e.BN.re.off <- sum((p.BN.re.off - test0$numerosi)^2)
e.ZI.Po.re.off <- sum((p.ZI.Po.re.off - test0$numerosi)^2)
e.ZI.BN.re.off <- sum((p.ZI.BN.re.off - test0$numerosi)^2)
errores_cuadraticos_re_AIC <- c(e.No.re, e.Po.re, e.Po.re.off, 
                                e.Po.re2, e.BN.re, e.ZI.Po.re, 
                                e.ZI.BN.re, e.Po.re2.off, e.BN.re.off, 
                                e.ZI.Po.re.off, e.ZI.BN.re.off)

Errores valor absoluto

Mostrar código
# Funcion de perdida valor absoluto
a.No.re <- sum(abs(p.No.re - test0$numerosi))
a.Po.re <- sum(abs(p.Po.re - test0$numerosi))
a.Po.re.off <- sum(abs(p.Po.re.off - test0$numerosi))
a.Po.re2 <- sum(abs(p.Po.re2 - test0$numerosi))
a.BN.re <- sum(abs(p.BN.re - test0$numerosi))
a.ZI.Po.re <- sum(abs(p.ZI.Po.re - test0$numerosi))
a.ZI.BN.re <- sum(abs(p.ZI.BN.re - test0$numerosi))
a.Po.re2.off <- sum(abs(p.Po.re2.off - test0$numerosi))
a.BN.re.off <- sum(abs(p.BN.re.off - test0$numerosi))
a.ZI.Po.re.off <- sum(abs(p.ZI.Po.re.off - test0$numerosi))
a.ZI.BN.re.off <- sum(abs(p.ZI.BN.re.off - test0$numerosi))
errores_abs_re_AIC <- c(a.No.re, a.Po.re, a.Po.re.off, a.ZI.BN.re, a.Po.re2, 
                        a.BN.re, a.ZI.Po.re, a.Po.re2.off,
                        a.BN.re.off, a.ZI.Po.re.off, a.ZI.BN.re.off)

Resumen errores (modelos con regularización) y AIC

Mostrar código
errores_re_AIC <- data.frame(
  modelo = c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", 
             "BN.re", "ZI.Po.re", "Po.re2.off",
             "BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off"),
  errores_cuadraticos_re_AIC,
  errores_abs_re_AIC)
errores_re_AIC[, 2:3] <- round(errores_re_AIC[, 2:3], 2)
errores_re_AIC %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()
modelo errores_cuadraticos_re_AIC errores_abs_re_AIC
No.re 270.56 192.77
Po.re 296.45 170.45
Po.re.off 458.35 203.40
ZI.BN.re 242.59 167.10
Po.re2 213.94 169.89
BN.re 222.83 167.08
ZI.Po.re 214.00 171.02
Po.re2.off 251.71 167.54
BN.re.off 213.78 163.79
ZI.Po.re.off 220.20 164.96
ZI.BN.re.off 213.78 163.80
Mostrar código
errores_combinados <- rbind(
  errores_re_AIC$errores_cuadraticos_re_AIC, 
  errores_re_AIC$errores_abs_re_AIC
)

colnames(errores_combinados) <- c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", "BN.re", "ZI.Po.re", "Po.re2.off",
             "BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off")

barplot(errores_combinados, 
        beside = TRUE, 
        names.arg = c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", 
             "BN.re", "ZI.Po.re", "Po.re2.off",
             "BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off"), 
        col = c(colores2[2], colores2[7]),
        border = "white",                
        main = "Comparación de Errores Cuadráticos y Absolutos", 
        xlab = "Modelo", 
        ylab = "Error", 
        las = 1,                      
        cex.names = 0.6,
        cex.axis = 0.8, 
        cex.main = 1.2)

Elegiríamos el modelo de la Binomial Negativa con offset. Entre los modelos evaluados, este modelo es el que presenta los mejores resultados generales, con el mínimo error tanto absoluto (163.79), como cuadrático (213.78). El valor del error cuadrático indica que el modelo maneja eficazmente posibles desviaciones extremas en los datos, mientras que el error absoluto refleja que en promedio sus predicciones se acercan más a los valores reales comparado con el resto de modelos.

Comparando con los modelos sin regularización, observamos que los errores de los modelos con regularización son más bajos en general. Existe una disminución de los errores en el modelo regularizado que hemos elegido respecto al modelo sin regularización anteriormente preferido.

4.2.4 Elección de modelos usando cross-validation

Mostrar código
train.m <- model.matrix(~., data=train)
train.m <- train.m[,colnames(train.m) != "numerosi"]

test.m <- model.matrix(~., data=test0)
test.m <- test.m[,colnames(test.m) != "numerosi"]
# Convertimos a matriz dispersa
train.m <- as(train.m, "dgCMatrix")

# Modelo de respuesta Normal
cv.No <- cv.glmnet(x=train.m, y=train$numerosi, alpha=0.5)

# Modelo de respuesta Poisson (glmnet)
cv.Po <- cv.glmnet(x=train.m, y=train$numerosi, alpha=0.5, 
                   family = "poisson")

# Modelo de respuesta Poisson y offset (glmnet)
cv.Po.off <- cv.glmnet(x = train.m, y = train$numerosi, 
                       alpha = 0.5, family = "poisson", 
                       offset = log(train$duracion_poliza/365.25))

# Modelo de respuesta Poisson (mpath)
#cv.Po2 <- cv.glmreg(numerosi ~ . , data = train, family="poisson",
#                    penalty = "enet")

# Modelo de respuesta Poisson y offset (mpath)

# Ajustar el modelo Poisson con penalización
#cv.Po2.off <- cv.glmreg(formula3 , data = train_log, 
#                        family = "poisson", penalty="enet",
#                        offset = log(duracion_poliza/365.25))

# Modelo de respuesta BN (mpath)
#cv.BN <- cv.glmregNB(formula3 , data= train, 
#                     penalty="enet", nfolds = 3, trace = TRUE)

# Modelo de respuesta BN y offset (mpath)
#cv.BN.off <- cv.glmregNB(formula3, data = train, 
#                         offset = log(train$duracion_poliza/365.25))

# Modelo de respuesta ZI-Po (mpath)
#cv.ZI.Po <- cv.zipath(numerosi ~ .| .  , data= train,
#                      family="poisson", link = "logit",
#                      penalty="enet")

# Modelo de respuesta ZI-BN (mpath)
#cv.ZI.BN <- cv.zipath(numerosi ~ .| . , data= train, 
#                      family="negbin", link = "logit",
#                      penalty="enet")

# Modelo de respuesta ZI-Po y offset (mpath)
# Variables y objetos auxiliares
train$duracion_poliza.a <- train$duracion_poliza/365.25
test0$duracion_poliza.a <- test0$duracion_poliza/365.25

formula5 <- formula(numerosi~ potencia +
                     antivehi+ tipovehi+ plazasvehi+ edad +
                     sexocondu+ anticia+
                     prov+ usovehi2 + g.potencia+
                     g.anticia+ edad2 + edad3 + edad_discretizada
                   +offset(log(duracion_poliza.a))
                   | potencia + antivehi+ tipovehi+
                     plazasvehi+ edad + sexocondu+
                     anticia+ prov+ usovehi2 + g.potencia+
                     g.anticia+ edad2 + edad3 + edad_discretizada
                   +offset(log(duracion_poliza.a))
)

#cv.ZI.Po.off<-cv.zipath(formula5, data= train, link = "logit",
#                        family= "poisson", penalty="enet")

# Modelo de respuesta ZI-BN y offset (mpath)
#cv.BN.Po.off <- cv.zipath(formula5, data= train, family= "negbin",
#                     link = "logit", penalty="enet")

El enfoque de cross-validation nos permite comparar el rendimiento de los modelos en diferentes conjuntos de entrenamiento y prueba, ayudando a identificar el modelo más adecuado según las métricas de error generadas durante la validación.

Predicciones con modelos ML con selección por CV

Mostrar código
p.cv.No <- predict(cv.No, test.m)

p.cv.Po <- predict(cv.Po, test.m, type="response")

p.cv.Po.off <- predict(cv.Po.off, test.m, type="response", 
                        newoffset = log(test$duracion_poliza/365.25))

#p.cv.Po2 <-predict(cv.Po2$fit, 
#     newx = test.m[, colnames(test.m) != "(Intercept)"], 
#     which= cv.Po2$lambda.which, type = "response")

#p.cv.Po2 <-predict(cv.Po2$fit, newx = test0)


#p.cv.Po2.off <- predict(cv.Po2$fit, 
#     newx = test.m[, colnames(test.m) != "(Intercept)"], 
#     which= cv.Po2$lambda.which, type = "response", 
#     newoffset = log(test0$duracion_poliza/365.25))

#p.cv.BN <-exp(predict(cv.BN, test0))

#p.cv.BN.off <- exp(predict(cv.BN.off, newx = test0, 
#                   newoffset = log(test0$duracion_poliza/365.25)))

#p.cv.ZI.Po<-predict(cv.ZI.Po, test0, 
#                    which= cv.ZI.Po$lambda.which, 
#                    type="response")

#p.cv.ZI.BN <-predict(ZI.BN.re, test0, 
#                     which= cv.ZI.BN$lambda.which, 
#                     type="response")

p.cv.ZI.Po.off<-predict(ZI.Po.off, newdata = test0, 
                        which=cv.ZI.Po.off$lambda.which,
                        type="response")

p.cv.ZI.BN.off<-predict(ZI.BN.off, newdata = test0, 
                        which=cv.ZI.BN.off$lambda.which,
                        type="response")

Tras entrenar los modelos utilizando técnicas de selección de hiperparámetros mediante validación cruzada, generamos las predicciones de los mismos.

Mostrar código
resultados <- data.frame(
  reales = test0$numerosi,
  cv_No = p.cv.No,
  cv_Po = p.cv.Po,
  cv_Po_off = p.cv.Po.off,
  cv_ZI_Po_off = p.cv.ZI.Po.off,
  cv_ZI_BN_off = p.cv.ZI.BN.off
)

resultados <- melt(resultados, id.vars = "reales", 
                   variable.name = "Modelo", 
                   value.name = "Predicciones")

resultados$Modelo <- recode(resultados$Modelo,
                            "lambda.1se" = "cv_No",
                            "lambda.1se.1" = "cv_Po",
                            "lambda.1se.2" = "cv_Po_off")

ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +
  geom_density(alpha = 0.4) +
  geom_density(data = resultados, aes(x = reales), 
               fill = "black", alpha = 0.2) +
  facet_wrap(~ Modelo) +
  scale_x_continuous(limits = c(0, 5)) +
  scale_fill_manual(values = colores3) +
  labs(x = "Número de Siniestros", 
       y = "Densidad", 
       title = "Comparación de Predicciones con Datos Reales") +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title.x = element_text(size = 12), 
    axis.title.y = element_text(size = 12),
    axis.text.x = element_text(),
    panel.border = element_blank(),
    plot.margin = margin(2, 1, 1, 1, "cm"),
    legend.position = "none")

Al igual que antes, el gráfico compara la densidad de las predicciones, esta vez de los modelos validados con validación cruzada. El modelo Binomial con offset muestra el mejor ajuste, logrando capturar de forma precisa la distribución de los datos reales, especialmente en los siniestros más frecuentes. En contraste, cv_No subestima los valores intermedios, y cv_Po_off tiende a sobreestimar siniestros altos, reduciendo su precisión. Los resultados refuerzan la preferencia por cv_ZI_BN_off, que demuestra una mejor capacidad de generalización.

Errores cuadráticos

Mostrar código
e.cv.No <- sum(abs(p.cv.No - test0$numerosi)^2)
e.cv.Po <- sum(abs(p.cv.Po - test0$numerosi)^2)
e.cv.Po.off <- sum(abs(p.cv.Po.off - test0$numerosi)^2)
#e.cv.Po2 <- sum(abs(p.cv.Po2 - test0$numerosi)^2)
#e.cv.Po2.off <- sum(abs(p.cv.Po2.off - test0$numerosi)^2)
#e.cv.BN <- sum(abs(p.cv.BN - test0$numerosi)^2)
#e.cv.BN.off <- sum(abs(p.cv.BN.off - test0$numerosi)^2)
#e.cv.ZI.Po <- sum(abs(p.cv.ZI.Po - test0$numerosi)^2)
e.cv.ZI.Po.off <- sum(abs(p.cv.ZI.Po.off - test0$numerosi)^2)
#e.cv.ZI.BN <- sum(abs(p.cv.ZI.BN - test0$numerosi)^2)
e.cv.ZI.BN.off <- sum(abs(p.cv.ZI.BN.off - test0$numerosi)^2)
errores_cuadraticos_re_CV <- c(e.cv.No, e.cv.Po, e.cv.Po.off,
                               e.cv.ZI.Po.off, e.cv.ZI.BN.off)

Errores valor absoluto

Mostrar código
# Funcion de perdida valor absoluto
a.cv.No <- sum(abs(p.cv.No - test0$numerosi))
a.cv.Po <- sum(abs(p.cv.Po - test0$numerosi))
a.cv.Po.off <- sum(abs(p.cv.Po.off - test0$numerosi))
#a.cv.Po2 <- sum(abs(p.cv.Po2 - test0$numerosi))
#a.cv.Po2.off <- sum(abs(p.cv.Po2.off - test0$numerosi))
#a.cv.BN <- sum(abs(p.cv.BN - test0$numerosi))
#a.cv.BN.off <- sum(abs(p.cv.BN.off - test0$numerosi))
#a.cv.ZI.Po <- sum(abs(p.cv.ZI.Po - test0$numerosi))
a.cv.ZI.Po.off <- sum(abs(p.cv.ZI.Po.off - test0$numerosi))
#a.cv.ZI.BN <- sum(abs(p.cv.ZI.BN - test0$numerosi))
a.cv.ZI.BN.off <- sum(abs(p.cv.ZI.BN.off - test0$numerosi))
errores_abs_re_CV <- c(a.cv.No, a.cv.Po, a.cv.Po.off, a.cv.ZI.Po.off, a.cv.ZI.BN.off)

Resumen errores (modelos con regularización) y CV

Mostrar código
errores_re_CV <- data.frame(
  modelo = c("No.re", "Po.re", "Po.re.off",
             "ZI.Po.re.off", "ZI.BN.re.off"),
  errores_cuadraticos_re_CV,
  errores_abs_re_CV)

errores_re_CV[, 2:3] <- round(errores_re_CV[, 2:3], 2)
errores_re_CV %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()
modelo errores_cuadraticos_re_CV errores_abs_re_CV
No.re 255.06 198.48
Po.re 234.75 184.17
Po.re.off 275.98 195.06
ZI.Po.re.off 382.59 192.51
ZI.BN.re.off 389.49 194.57
errores_combinados <- rbind(
  errores_re_CV$errores_cuadraticos_re_CV, 
  errores_re_CV$errores_abs_re_CV
)

colnames(errores_combinados) <- c("No.re", "Po.re", "Po.re.off","ZI.Po.re.off", "ZI.BN.re.off")
Mostrar código
barplot(errores_combinados, 
        beside = TRUE, 
        names.arg = c("No.re", "Po.re", "Po.re.off","ZI.Po.re.off", "ZI.BN.re.off"), 
        col = c(colores2[2], colores2[7]),
        border = "white",                
        main = "Comparación de Errores Cuadráticos y Absolutos", 
        xlab = "Modelo", 
        ylab = "Error", 
        las = 1,                      
        cex.names = 1,
        cex.axis = 0.8, 
        cex.main = 1.2)

Resumen errores modelos con regularización

Mostrar código
errores_re <- cbind(errores_re_AIC[c(1:3,9,8),], errores_re_CV[,-1])
errores_re %>%
  kbl(align = rep('c', 3)) %>%
  kable_styling()
modelo errores_cuadraticos_re_AIC errores_abs_re_AIC errores_cuadraticos_re_CV errores_abs_re_CV
1 No.re 270.56 192.77 255.06 198.48
2 Po.re 296.45 170.45 234.75 184.17
3 Po.re.off 458.35 203.40 275.98 195.06
9 BN.re.off 213.78 163.79 382.59 192.51
8 Po.re2.off 251.71 167.54 389.49 194.57

Volvemos a elegir el modelo BN.re.off como la mejor opción. Aunque fijándonos únicamente en el valor de los errores el modelo No.re podría ser más adecuado, este no incluye mecanismos para manejar sobredispersión ni incorpora términos de regularización que ayuden a evitar el sobreajuste

El modelo de la Binomial Negativa es particularmente adecuado para nuestro conjunto de datos, ya que, como se ha mencionado, maneja eficientemente la sobredispersión que caracteriza a la variable dependiente. Además, la regularización, a través del término offset, contribuye a evitar problemas de sobreajuste y mejora la estabilidad del modelo.

Comparando con los modelos anteriores, BN.re.off demuestra una mejora significativa respecto a los modelos sin regularización, como BN.log.off, que aunque presenta buenos resultados, no alcanza el desempeño del modelo seleccionado en las métricas clave. En relación a los modelos evaluados únicamente con AIC (evaluación clásica), BN.re.off logra un balance superior entre ajuste y capacidad predictiva, superando a opciones como ZI.Po.off y ZI.BN, cuyos valores de AIC más bajos no se traducen en un mejor desempeño en validación cruzada.

5. Conclusiones

En conclusión, dado que nuestro estudio se basa en una muestra donde la varianza de los datos supera la media, la Binomial Negativa es preferible gracias a su capacidad para capturar la estructura de estos datos, manejar la sobredispersión y garantizar predicciones confiables en nuevos conjuntos.

La decisión de reducir la muestra nos permitió agilizar la convergencia de los modelos y continuar con el análisis. Sin embargo, es importante destacar que esto podría haber limitado la capacidad de los modelos para captar patrones más complejos presentes en la muestra completa. Adicionalmente, se eliminaron algunas observaciones específicas para evitar problemas al realizar las predicciones, como aquellas correspondientes a niveles de factores no presentes en el conjunto de entrenamiento. Estos ajustes fueron necesarios para garantizar la ejecución, pero podrían haber afectado la representatividad de los resultados. Por lo tanto, sería interesante repetir el análisis utilizando la totalidad de los datos originales en un equipo de soporte más preparado para manejar tal cantidad de datos.

6. Referencias

[Callout Blocks – quarto. (n.d.). Quarto.] (https://quarto.org/docs/authoring/callouts.html)

[INE - Instituto Nacional de Estadística. (n.d.). INEbase/ Clasificaciones /Relación de municipios, provincias, comunidades y ciudades autónomas con sus códigos / Relación de comunidades y ciudades autónomas con sus códigos. INE - Instituto Nacional De Estadística.] (https://www.ine.es/daco/daco42/codmun/cod_ccaa_provincia.htm)

[¿Qué es la correlación estadística y cómo interpretarla? (n.d.). Máxima Formación.] (https://www.maximaformacion.es/blog-dat/que-es-la-correlacion-estadistica-y-como-interpretarla/)

[colaboradores de Wikipedia. (2023, November 22). Distribución log-normal. Wikipedia, La Enciclopedia Libre.] (https://es.wikipedia.org/wiki/Distribuci%C3%B3n_log-normal)

[Color palettes – Data Visualization with R. (n.d.). Data Visualization With R.] (https://datavizs24.classes.andrewheiss.com/resource/colors.html)

[What are continuous probability distributions & their 8 common types? | KNIME. (n.d.). KNIME.] (https://www.knime.com/blog/continuous-probability-distribution)